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I.  INTRODUCTION 


A.  BACKGROUND 

The  development  of  rotorcraft  design  is  a  complicated  iterative  process  involving 
the  interaction  of  many  variables.  It  is  an  especially  difficult  and  time  consuming 
task  when  the  required  calculations  are  performed  by  hand.  There  are  several 
helicopter  design  software  packages  that  perform  the  required  calculations  and 
automate  the  design  analysis  process.  In  searching  for  a  computer  program  suitable 
for  helicopter  preliminary  design,  students  found  that  current  helicopter  design  codes 
were  either  too  difficult  to  use,  excessively  detailed,  oi  iiuiwuonally  incomplete.  As  a 
result,  the  need  for  a  custom  computer  code  was  identified.  Through  this  thesis,  the 
Joint  Army/Navy  Rotorcraft  Analysis  and  Design  (JANRAD)  computer  program  was 
developed  to  meet  the  specific  needs  of  the  Naval  Postgraduate  School  for  preliminary 
helicopter  design. 

JANRAD  was  written  as  an  interactive,  user  friendly  program,  capable  of 
accurately  and  quickly  solving  helicopter  design  problems  at  the  preliminary  design 
level.  The  program  was  intended  primarily  for  analysis  of  helicopter  rotor 
performance,  stability  and  control,  and  rotor  dynamics.  It  was  not  intended  to 
perform  comprehensive  helicopter  design  analyses. 

JANRAD  was  written  as  a  collection  of  MATLAB*  script  files  (M-files).  It  was 
written  specifically  for  version  3,5  of  386-MATLAB*  [Ref,  1].  This  was  done  for 
two  reasons.  First,  MATLAB*  was  ideally  suited  for  the  types  of  calculations  and 
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analyses  involved  in  the  preliminary  helicopter  design  process.  Second,  engineering 
design  students  had  typically  used  MATLAB*  for  several  projects  and  had  a  good 
working  knowledge  of  the  software  and  the  programming  language.  The  intent  was  to 
have  a  powerful  and  accurate  computer  aided  design  tool  that  one  could  master  in  a 
matter  of  minutes.  In  using  JANRAD,  students  could  occupy  their  time  analyzing 
design  issues  rather  than  learning  software. 

B.  APPROACH  TO  THE  PROBLEM 

The  approach  in  developing  JANRAD  was  to  provide  the  design  student  with  a 
rapid  means  of  changing  design  parameters  and  obtaining  accurate  results  over  a  wide 
range  of  helicopter  configurations.  To  accomplish  this,  user  input  requirements  were 
kept  to  a  minimum,  efficiency  in  parameter  calculation  routines  was  optimized, 
simplifying  assumptions  were  made,  and  convergence  criteria  for  iterative  processes 
were  carefully  set  to  appropriate  levels. 

To  keep  user  input  to  a  minimum,  airframe  and  rotor  geometric  input 
parameters  were  minimized.  The  rotor  blade  was  limited  to  a  rectangular,  uniform, 
rigid  structure  with  linear  twist.  This  allowed  single  value  inputs  for  setup  of  the 
rotor  blade  parameters.  The  wing  and  tail  section  inputs  were  limited  by  using  area 
and  span  to  determine  other  geometric  values.  User  input  requirements  were  also 
reduced  by  requiring  aircraft  gross  weight,  equivalent  flat  plate  area,  auxiliary  thrust, 
and  aerodynamic  surface  coefficients  of  lift  to  be  precalculated. 

To  maximize  the  efficiency  of  the  iterative  processes,  whenever  possible, 
MATLAB’ s*  built-in  vector  and  matrix  operations  were  used  in  place  of  repetitive 
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"for"  and  "while"  loc^s.  In  addition,  adjustments  to  the  indq)endent  variables  during 
the  iterative  processes  were  varied  according  to  certain  parameters  such  as  advance 
ratio,  density  ratio,  and  magnitude  of  error.  By  varying  the  increment  or  decrement 
of  the  independent  variables  according  to  the  appropriate  parameters,  the  speed  of 
convergence  was  greatly  increased. 

In  making  assumptions  and  setting  convergence  criteria,  the  trade-offs  between 
speed  and  accuracy  were  balanced  by  keeping  in  mind  that  JANRAD  was  developed 
for  helicopter  preliminary  design.  Accuracy  was  maintained  for  a  wide  range  of 
helicopter  configurations  by  using  blade  element  theory  to  determine  resultant  forces 
and  moments  in  the  rotor  system.  Overall  speed  was  maintained  by  making 
simplifying  assumptions  and  selecting  convergence  criteria  that  resulted  in  small 
percentage  errors  that  were  acceptable  for  preliminary  design. 

C.  PROGRAM  OUTLINE 

Currently,  JANRAD  consists  of  a  group  of  22  MATLAB*  script  and  function 
files.  The  main  IA-T\\c,  janrad.m,  invokes  the  user  interface  functions  and  the  various 
analysis  modules.  The  routine  allows  the  user  to  interactively  input  environment  and 
design  configuration  parameters  using  the  computer  keyboard  and  screen.  The  user 
can  either  create  or  edit  files  containing  parameters  for  specific  designs  and  situations. 
The  routine  also  provides  the  interface  for  selection  of  the  various  design  analyses.  A 
listing  of  janrad.m  is  found  at  Appendix  A.  Currently,  there  is  a  module  available 
for  rotor  performance  and  one  for  stability  and  control.  A  rotor  dynamics  routine  is 
under  development. 
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The  heart  of  JANRAD  is  the  M-file  trim.m.  This  M-fUe  receives  the 
environment  and  (tesign  configuration  parameters  from  janrad.m  and  then  balances  the 
forces  and  moments  in  the  rotor  system  for  the  appropriate  airspeed.  By  trimming 
the  forces  in  the  rotor  system,  trim.m  produces  the  parameters  necessary  for  the 
various  analysis  modules.  A  listing  of  trim.m  is  found  at  A|^ndix  B. 

The  rotor  performance  module  is  called  perf.m.  This  M-file  accepts  the  data 
from  trim.m  and  janrad.m,  calculates  the  performance  parameters,  and  generates  the 
rotor  performance  output.  Appendix  C  contains  a  listing  of  perf.m.  The  stability  and 
control  module,  stab.m  calculates  and  produces  the  stability  and  control  output.  Each 
M-file,  with  the  exception  of  stab.m,  is  discussed  in  detail  in  this  thesis.  The  details 
of  the  stability  and  control  module  are  addressed  in  another  thesis  [Ref.  2]. 

Similarly,  the  details  of  the  rotor  dynamics  module  will  be  discussed  in  a  future 
thesis. 
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n.  PROGRAM  DETAILS 


A.  GENERAL 

JANRAD  provides  accurate  data  for  standard  or  compound  helicopters  in 
hovering  and  forward  flight.  Forward  flight  analysis  is  only  accurate  for  airspeeds  of 
50  knots  or  greater.  Conditions  are  limited  to  level  flight  with  no  winds.  JANRAD 
provides  data  for  one  flight  condition  and  one  design  configuration  at  a  time.  To 
obtain  a  range  of  data  for  an  overall  analysis,  the  user  must  make  multiple  runs  of 
JANRAD.  For  each  run,  the  appropriate  parameters  would  be  changed  and  the 
results  could  then  be  compared  to  those  of  previous  runs. 

B,  THEORY 

The  underlying  principle  in  helicopter  rotor  design  is  that  the  main  rotor  forces 
and  moments  must  be  adjusted  so  the  rotor  provides  the  lift  and  propulsive  thrust 
required  to  meet  the  conditions  of  flight.  The  process  of  adjusting  the  rotor’s  forces 
and  moments  is  referred  to  as  "trimming"  the  rotor.  At  a  hover  with  no  wind  and 
assuming  centered  C.G.  and  no  rolling  moment  due  to  the  anti-torque  device,  there  is 
no  tilt  required  in  the  rotor  disk.  In  this  case,  trimming  the  rotor  involves  adjusting 
collective  pitch  until  the  thrust  generated  by  the  rotor  equals  the  weight  of  the  aircraft 
plus  the  force  created  by  the  rotor  downwash  on  the  fuselage.  In  level,  unaccelerated 
forward  flight,  assuming  a  centered  C.G.  and  no  rolling  moments,  the  rotor  disk  must 
be  tilted  forward  to  create  a  vertical  and  a  horizontal  thrust  component.  The  vertical 
thrust  component  must  equal  the  weight  of  the  aircraft  (no  downwash  effects)  minus 
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lift  created  by  any  aerodynamic  surfaces  (assuming  no  lift  generated  by  the  fuselage). 
The  horizontal  component  must  equal  the  drag  created  by  the  fuselage  and  other 
aircraft  components  minus  any  auxiliary  thrust.  Trimming  the  rotor  in  forward  flight 
involves  the  adjustment  of  cyclic  pitch  (collective  pitch  is  the  steady  term  of  cyclic 
pitch)  to  maintain  the  required  rotor  disk  tilt. 

For  hovering  and  forward  flight,  JANRAD  uses  blade  element  theory  to 
determine  the  rotor  system  forces  and  moments,  and  a  method  for  harmonically 
balancing  the  forces  and  moments  to  adjust  cyclic  pitch  [Ref.  3].  The  only  difference 
in  the  process  for  trimming  the  rotor  system  in  hovering  and  forward  flight  is  in  how 
the  rotor  induced  velocity  is  calculated.  In  forward  flight,  the  rotor  induced  velocity 
is  assumed  uniform  and  is  calculated  using  the  inflow  parameter  Xj  [Ref.  4].  For 
hovering  flight,  induced  velocity  is  determined  by  using  an  iterative  process  of 
equating  the  differential  thrust  calculated  by  momentum  theory  and  blade  element 
theory  [Ref.  5,  p.  38].  This  process  produces  an  elliptical  non-uniform  induced 
velocity  distribution  along  the  rotor  blade. 

C.  ASSUMPTIONS 

The  following  assumptions  are  made  in  order  to  reduce  total  computation  time. 
The  assumptions  produce  relatively  small  errors  which  are  considered  acceptable  for 
helicopter  preliminary  design. 

1.  Center  of  Gravity 

For  the  purposes  of  trimming  the  rotor  system,  the  aircraft  center  of 
gravity  is  located  at  the  center  of  the  rotor  hub.  In  doing  so,  the  first-harmonic  terms 
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of  thrust  moment  are  eliminated. 

2.  Lateral  Forces  and  Moments 

The  rotor  is  trimmed  longitudinally  and  vertically.  Lateral  forces  and 
moments  due  to  the  anti-torque  device  and  aerodynamic  surfaces  are  ignored. 

3.  Rigid  Rotor  Blade 

For  the  purposes  of  the  rotor  trimming  routine,  the  rotor  blade  is 
considered  a  rigid  structure.  This  eliminates  the  need  to  consider  blade  modal 
responses  and  the  subsequent  affects  on  the  rotor  blade  local  angle  of  attack. 

4.  Tip  Loss 

Rotor  blade  tip  loss  is  determined  using  a  relationship  involving  coefficient 
of  thrust.  This  results  in  a  steady  value  of  tip  loss  with  respect  to  azimuth  angle. 

The  harmonic  characteristics  of  tip  loss  are  not  considered. 

5.  Aerodynamic  Surfaces 

Lift  and  drag  for  the  aircraft  wing,  horizontal  tail,  and  vertical  tail  are 
calculated  using  two  dimensional  theory.  Spanwise  flow  and  tip  losses  are  not 
accounted  for.  The  aircraft  fuselage  lift  characteristics  are  also  ignored. 

6.  Mach  Number 

The  airfoil  lift  and  drag  curves  (coefficients  of  lift  and  drag  versus  angle  of 
attack)  for  separate  Mach  numbers  were  reduced  to  one  Mach  independent  curve. 

7.  Wake  Interactions 

Variations  in  angle  of  attack  due  to  the  interaction  of  rotor  wake, 
aerodynamic  surface  trailing  vortices,  or  fuselage  interference  were  not  accounted  for. 
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D.  M-FILE  DETAILS 


1.  Janrad.m 

The  main  module  for  the  JANRAD  program  is  janrad.m.  This  M-file 
first  queries  the  user  to  either  edit  or  create  a  data  file.  If  editing  is  selected,  the 
program  asks  for  the  file  name.  If  the  file  name  is  not  found  in  the  current  directory, 
the  program  will  tell  the  user  to  try  again  or  instruct  them  on  how  to  exit  JANRAD. 
When  the  file  name  is  found,  a  list  of  parameters  is  displayed,  and  the  user  is  allowed 
to  make  changes  as  desired.  If  the  user  decides  to  create  a  new  data  file,  the 
parameters  are  displayed  one  at  a  t'*^e.  In  each  instance,  the  user  is  prompted  to 
input  the  desired  value.  Each  prompt  includes  a  display  of  the  proper  units.  The 
parameters  eligible  for  editing  or  creation  are  identical  and  are  listed  below: 

1 .  pressure  altitude 

2.  temperature 

3.  airspeed 

4.  aircraft  gross  weight 

5.  number  of  blades  in  the  rotor  system 

6.  rotor  blade  radius 

7.  rotor  blade  chord 

8.  effective  hinge  offset 

9.  rotor  blade  grip  length 

10.  rotor  blade  angle  of  twist 

1 1 .  rotor  blade  weight 

12.  number  of  rotor  blade  elements 

13.  rotor  system  rotational  velocity 

14.  number  of  rotor  disk  azimuth  sectors 

15.  average  lift  curve  slope  of  rotor  blade  airfoil 

16.  rotor  blade  airfoil 

17.  collective  pitch  angle  at  70%  rotor  blade  radius 

18.  equivalent  aircraft  flat  plate  area 

19.  aircraft  vertical  projected  area 

20.  wing  area 

21.  wing  span 

22.  wing  lift  coefficient 
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23.  wing  proflle  drag  coefficient 

24.  wing  efficiency  factor 

25.  horizontal  tail  area 

26.  horizontal  tail  span 

27.  horizontal  tail  lift  coefficient 

28.  horizontal  tail  profile  drag  coefficient 

29.  vertical  tail  area 

30.  vertical  tail  span 

31.  vertical  tail  lift  coefficient 

32.  vertical  tail  profile  drag  coefficient 

33.  auxiliary  thrust 

Once  the  editing  or  creation  process  is  complete,  the  program  then  queries 
the  user  for  a  file  name.  JANRAD  saves  the  data  to  the  current  directory  with  the 
specified  file  name  and  a  ".mat"  extension.  This  is  a  binary  file  containing  the  file 
label  and  the  names  and  values  of  the  33  parameters. 

Finally,  the  user  is  queried  for  a  specific  analysis.  To  date,  the  user  can 
select  rotor  performance  or  aircraft  stability  and  control.  Janrad  has  provisions  for  a 
rotor  system  dynamics  module.  When  the  user  selects  rotor  performance,  the  module 
peif.m  is  called.  Similarly,  for  stability  and  control,  stab.m  is  called.  The  user  can 
also  make  selections  to  either  change  the  current  data  (edit  current  parameters  or 
create  new  ones)  or  exit  the  program. 

2.  Trim.m 

The  M-file  trim.m  trims  the  rotor  system  as  appropriate  for  the 
environment  and  design  parameters  obtained  as  input  in  janrad.m.  The  rotor 
trimming  process  involves  four  steps: 

1 .  Calculate  required  parameters. 

2.  Guess  rotor  drag,  cyclic  pitch,  and  location  of  resultant  thrust  vector. 

3.  Trim  the  rotor  system  (calculate  cyclic  pitch). 

4.  Calculate  rotor  drag  and  location  of  resultant  thrust  vector. 
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Initially,  the  prognun  uses  values  for  rotor  drag,  cyclic  pitch,  and  the  location  of  the 
resultant  thrust  vector  based  on  established  averages.  The  rotor  system  is  then 
trimmed  using  those  values  and  the  other  parameters  previously  calculated.  Once  the 
rotor  is  trimmed,  rotor  drag  and  the  location  of  the  resultant  thrust  vector  are 
calculated.  If  the  values  for  these  two  parameters  fall  within  certain  evaluation 
criteria  (20%  of  prior  rotor  drag  and  1.5%  of  prior  location  of  resultant  thrust 
vector),  then  the  rotor  is  considered  trimmed.  If  not,  the  rotor  is  retrimmed  using 
new  values  for  rotor  drag  and  the  location  of  the  resultant  thrust  vector.  This  process 
is  repeated  until  the  two  values  fall  within  the  evaluation  criteria. 

As  the  first  step,  trim.m  calculates  various  parameters  used  throughout  the 
routine.  The  parameters  include  ambient  air  density  (pj,  rotor  drag  or  H-force 
(Hr«of)>  dynamic  pressure  (q),  rotor  disk  area  rotor  blade  tip  speed  (V„p), 
fuselage  drag  lift  and  drag  on  the  aerodynamic  surfaces  (D,i„g,  I^veit 

<iia>  Lwing.  Lhorizuu*  Uenuii).  the  lotor  tip  path  plane  angle  (otpp),  and  the  advance  ratio 
(ji).  The  calculation  of  rotor  drag  at  this  point  is  a  first  guess  based  on  the 
relationship; 

(1) 

y  o 

where: 

V.  =  forward  airspeed 
p,  =  ambient  air  density 
Po  =  sea  level  density 
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The  calculations  for  lift  and  drag  on  the  wing  and  tail  sections  are  based  on  a  two 
dimensional  analysis  where: 


CL 

=  c„  + 

(2) 

It  eAR 

D  —  (jCpS 

(3) 

L  =  qC^S 

(4) 

where: 

Cd  =  drag  coefficient 
Cl  =  lift  coefficient 
Cdo  =  profile  drag  coefficient 
e  =  wing  efficiency  factor 
AR  =  wing  aspect  ratio 
q  =  dynamic  pressure 
S  =  wing  area 

The  values  for  the  horizi.  r.tal  and  vertical  tail  efficiency  factors  are  set  at  0.8.  The 
efficiency  factor  for  the  wing,  as  well  as  Cl,  and  Cqo  for  the  wing  and  other  surfaces 
are  direct  input  parameters  obtained  from  janrcui  m.  The  drag  coefficients,  aspect 
ratios,  and  areas  of  the  aerodynamic  surfaces,  as  well  as  q  are  calculated  from  inpnt 
parameters  obtained  from  janrad.m.  The  rotor  tip  path  plan  angle  is  calculated  as: 
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(5) 


Oj.pp 


=  tan"^ 


GW-Lt  ) 


where: 

Df  —  +  D^ipf  ■t"  D|io^  +  ^vcft  uti '  T,iu 

=  auxiliary  thrust 
Hmot  =  H-force  (rotor  drag) 

GW  =  aircraft  gross  weight 

Lf  Lwii^  "i" 

Finally,  the  advance  ratio  is  calculated  as: 

_  U.COS  (ttypp) 


(6) 


where: 

V.  =  forward  velocity 
otpp  =  tip  path  plane  angle 
Vyp  =  rotor  blade  tip  speed 

The  required  thrust  and  corresponding  coefficient  of  thrust  also 
calculated.  For  hovering  flight,  thrust  required  is  calculated  as: 


r= 


1+0.3 


Voxt  pzoj 
^disk  , 


GF/ 


where: 

Avert  pK!j  =  vertical  projected  area 
Ajirt  =  area  of  rotor  disk 


(7) 
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GW  =  aircraft  gross  weight 

This  equation  assumes  an  effective  vertical  drag  coefficient  of  0.3  for  all  the  aircraft 
components  in  the  rotor  wake  [Ref.  5,  p.  7].  In  forward  flight,  the  thrust  required  is 
calculated  as: 


GW-Lf 
COS  (ffypp) 


(8) 


where: 


GW  =  aircraft  gross  weight 

f^ing  f'horiz  Uil 

a^pp  =  tip  path  plane  angle 
The  coefficient  of  thrust  is  found  by: 


Cj.= 


^disk  Pa  ^tip 


(9) 


where: 


T  =  thrust 

Adi*  =  rotor  disk  area 
p,  =  ambient  air  density 
V,ip  =  rotor  blade  tip  speed 

Once  Cf  is  calculated,  the  rotor  blade  tip  loss  factor  is  determined  as: 


B  =  l- 


2C, 


(10) 


13 


where: 

Ct  =  coefficient  of  thrust 

b  =  number  of  blades 

Once  the  required  parameters  are  calculated,  trim.m  sets  up  the  rotor  blade 
radial  sections  and  rotor  disk  azimuth  sectors.  The  radial  sections  are  determined  by 
evenly  dividing  the  distance  from  the  beginning  of  the  aerodynamic  portion  of  the 
rotor  blade  to  the  effective  blade  radius  (Refr  =  B'*'R)  by  the  number  of  blade  elements 
input  by  the  user  in  janrad.m.  Azimuth  is  determined  by  dividing  the  360  degree 
rotor  disk  into  equally  spaced  areas  determined  by  the  number  of  azimuth  sectors 
input  by  the  user. 

The  program  then  calculates  rotor  blade  geometric  angle  due  to  twist  (jSJ 
and  rotor  coning  angle  Blade  twist  is  considered  linear  and  determined  by  the 
relationship: 

Pt=e.(o.7-^)  (11) 

Where: 

=  blade  twist 

r  =  distance  from  center  of  hub  to  center  of  blade  element 

R  =  rotor  blade  radius 

Equation  1 1  takes  this  form  because  blade  twist  is  converted  to  a  positive  value  when 
input  by  the  user  in  janrad.m.  This  allows  the  user  to  input  either  a  positive  or 
negative  value  of  blade  twist  without  affecting  the  accuracy  of  the  program’s  output. 
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It  would  be  highly  unlikely  that  blade  twist  would  ever  have  positive  value  yet  some 
may  refer  to  it  as  positive  when  it  is  actually  negative.  The  rotor  coning  angle  is 
d^rmined  by: 


Po  =  sin-^ 


(12) 


where: 

$i  =  blade  twist 

r  =  distance  from  center  of  hub  to  center  of  blade  element 
T  =  thrust 

b  =  number  of  blades 
R?  =  (R,ff  -  e)’''rT 

R  =  rotor  blade  radius 
e  =  effective  hinge  offset 
Tt  =  location  of  resultant  thrust  vector 
R,ff  =  effective  radius  (B*R) 

Wbude  =  weight  of  rotor  blade 
0  =  rotor  rotational  velocity 
nibude  =  mass  of  rotor  blade 

Equation  12  is  based  on  a  uniform  rectangular  blade,  uses  small  angle  approximations 
for  jSo,  and  initially,  assumes  the  location  of  the  resultant  thrust  vector  is  at  0.7  r/R. 
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The  next  stq>  is  to  determine  the  induced  velocity  distribution  (vj.  At  a 
hover,  Vi  is  determined  by  equating  differential  thrust  from  momentum  theory: 

dr=4  r  dr  (13) 

and  that  of  blade  element  theory: 

dr*jb-^  (Qr)2a^e--^jcdr  d*) 

which  results  in  a  quadratic  equation  of  the  form: 

bac^v 

where: 

p,  =  ambient  air  density 
0  =  rotor  rotational  velocity 
b  =  number  of  blades 
a  =  airfoil  lift  curve  slope 
c  =  rotor  blade  chord 
6  =  blade  pitch  angle  (/3,  +  6^) 

/3,  =  geometric  angle  due  to  twist 
0^  =  collective  pitch  angle 

r  =  distance  from  center  of  hub  to  center  of  blade  element 
dr  =  blade  element  width 

The  maximum  root  for  Vj  is  then  taken  and  used  in  equation  14  to  determine  the 
differential  thrust.  These  calculations  are  repeated,  incrementing  $  appropriately  until 


(47t)  vj 


^-i  —  rbacB  *0 
M2 


(15) 
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the  sum  of  the  dT’s  equals  the  required  thrust  determined  from  equation  7.  In 
forward  flight,  induced  velocity  is  considered  uniform  and  is  calculated  using  the 
inflow  parameter  (X^)  where: 


Xj.=  p  sin(a„p)  + 


—  C 


A 


The  maximum  root  for  X^  is  taken  and  Vi  is  calculated  as: 


=Xj.Qi?-V.sin(aj.pp) 


(16) 


(17) 


where: 

ft  =  advance  ratio 
qtpp  =  tip  path  plane  angle 
Cj  =  coefficient  of  thrust 
Q  =  rotor  rotational  velocity 
R  =  rotor  blade  radius 
V»  =  forward  airspeed 

The  last  step  before  starting  the  rotor  trim  routine  is  to  make  a  first  guess 
at  cyclic  pitch  (6).  Cyclic  pitch  is  a  first  harmonic  function  with  the  cosine  term  (©iJ 
having  a  small  positive  value  dependent  on  airspeed  and  the  sine  term  (^,,),  a  negative 
value,  also  dependent  on  airspeed.  The  steady  term  or  collective  pitch  term  is  an 
input  parameter  obtained  from  the  user  in  janrad.m.  The  cyclic  pitch  is  calculated  as: 

0=0o+0ieCos  (ijr) +0i5Sin(i)r)  (18) 


where:  \J/  =  azimuth  angle 


The  rotor  system  is  trimmed  using  blade  element  theory  in  a  three  step 
iterative  process  that  consists  of  first  adjusting  collective  pitch  (dj,  then  trimming 
cyclic  pitch  (0),  and  finally  retrimming  the  collective.  Blade  element  theory  is  based 
on  the  determination  of  local  angle  of  attack  (a)  for  each  rotor  blade  element.  Once 
a  is  calculated,  the  coefficient  of  lift  (C  J  and  coefficient  of  drag  (Cq)  for  each 
element  are  determined  and  fmally,  thrust  and  drag  for  each  element  are  calculated. 
Blade  element  theory  involves  the  use  of  the  following  equations: 

o  =  0+p£-«l)  (19) 


<j>  =  tan‘^ 


(20) 


C7p  =  v^cos(Po)  +v:sin(a„p)cos(po)  .  . 

+  ^.008  (ffypp)  sin(pj,)  cos  (ilr) 

C/j.  =  Qr  +  V.cos  (Oj-pp)  sin(ilr)  (22) 

where: 

$  =  cyclic  pitch  (eqn.  18) 

P,  =  rotor  blade  geometric  angle 
<f>  =  inflow  angle 

Up  =  velocity  perpendicular  to  rotor  blade 
U,  =  velocity  parallel  to  rotor  blade 
Vj  =  induced  velocity 
Po  =  coning  angle 
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V»  =  forward  airspeed 
aypp  =  tip  path  plane  angle 
^  =  azimuth  angle 

r  =  distance  from  center  of  hub  to  blade  element 
The  rotor  trim  procedure  use  two  M-files,  thrcalc.m  and  rmcalc.m.  The 
listings  for  these  M-files  are  found  at  Appendicies  D  and  E  respectively.  The  first  M- 
file,  thrcalc.m,  uses  blade  element  theory  to  calculate  the  differential  thrust  along  the 
rotor  blade  at  each  azimuth  angle.  The  differential  thrust  along  the  rotor  blade  is  then 
summed  to  give  the  total  thrust  at  each  azimuth  angle  (T^).  Differential  thrust  is 
calculated  as: 

dT^^^p^cdr{ul*ul)  [qcos(<|>)-C^sin(<|))]  (23) 

where; 

p,  =  ambient  air  density 
c  =  rotor  blade  chord 
dr  =  blade  element  width 
Up  =  velocity  perpendicular  to  rotor  blade 
U,  =  velocity  parallel  to  rotor  blade 
Cl  =  coefficient  of  lift 
Cd  =  coefficient  of  drag 
=  inflow  angle 

The  small  contribution  of  thrust  provided  by  profile  drag  in  the  tip  loss  region  of  the 


19 


blade  is  accounted  for  by  giving  the  profile  drag  coefficient  (CdJ  a  nominal  value  of 
.009.  The  equation  for  differential  thrust  in  the  tip  loss  region  is  as  follows; 

I  P.C  ( - .  009  sin  (♦„,) )  (24) 

where: 

p,  =  ambient  air  density 

c  =  rotor  blade  chord 

R  =  rotor  blade  radius 

R,ff  =  effective  rotor  radius 

Up  pp  =  velocity  peipendicular  to  rotor  blade 

U,iip  =  velocity  parallel  to  rotor  blade 

0  =  inflow  angle 

The  second  M-file,  tmcalc.m,  also  uses  blade  element  theory  to  calculate  the 
differential  thrust  moments  along  the  blade  at  each  azimuth  angle.  The  differential 
thrust  moments  are  then  summed  to  give  the  total  thrust  moment  at  each  azimuth 
angle  (M^).  The  differential  thrust  moment  is  calculated  as: 

dM^=^p^crdr[ul+ul)  [C^cos  (4>)  -CpSin((|>)  ]  (25) 

where: 

p,  =  ambient  air  density 
c  =  rotor  blade  chord 

r  =  distance  from  center  of  hub  to  blade  element 
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dr  =  blade  element  width 


Up  =  velocity  perpendicular  to  rotor  blade 
U,  =  velocity  parallel  to  rotor  blade 
Cl  =  coefficient  of  lift 
Co  =  coefficient  of  drag 
<f>  =  inflow  angle 

The  small  thrust  moment  provided  by  profile  drag  in  the  tip  loss  region  of  the  blade  is 
accounted  for  in  a  similar  manner  as  the  thrust  in  thrcalc.m.  The  equation  is: 

1  (  (R-R  )  ^  (26) 

dMtix>=\  ( - .  009  sin  ) 

where: 

p,  =  ambient  air  density 

c  =  rotor  blade  chord 

R  =  rotor  blade  radius 

R,ff  =  effective  rotor  radius 

Up  Up  =  velocity  perpendicular  to  rotor  blade 

U,up  =  velocity  parallel  to  rotor  blade 

<f>  =  inflow  angle 

The  two  M-files,  thrcalc.m  and  tmcalc.m,  caU  an  M-file  which  contains  the 
lift  coefficient  and  the  drag  coefficient  data  for  the  appropriate  rotor  blade  airfoil. 

The  airfoil  data  are  arranged  as  separate  polynomials  for  various  ranges  of  angle  of 
attack  (a).  There  are  currently  two  airfoil  data  function  files,  hh02clcd.m  and 


wl2clcd.m.  The  former  contains  data  for  the  McDonnell  Douglas  Helio^ter 
Corporation  (MDHC)  HH-02  airfoil  used  on  the  AH-64A  Apache  attack  helicopter. 
The  latter  contains  data  for  the  Boeing  VR-12  airfoil  planned  for  use  on  the  RAH-66 
Comanche.  There  is  no  Mach  number  dq)endency  incorporated  in  either  airfoil  data 
function  flies.  The  listings  for  these  M-files  are  found  in  Appendicies  F  and  G 
respectively. 

In  both  the  initial  trim  and  retrim  portions  of  the  coUective  trim  procedure, 
the  collective  is  adjusted  by  increasing  or  decreasing  6^  a  suitable  amount  to  match  the 
required  thrust  calculated  from  either  equation  7  or  8.  The  cyclic  terms  (0,c  and 
are  not  adjusted  during  the  collective  trim  procedure.  Collective  pitch  is  adjusted 
until  the  resulting  thrust  calculated  by  thrcalc.m  is  within  2%  of  required  thrust  for 
the  initial  trim  procedure  and  1  %  of  required  thrust  for  the  retrim  procedure. 

The  cyclic  trim  routine  begins  by  determining  the  initial  thrust  moment 
based  on  the  cyclic  pitch  determined  by  the  collective  trim  procedure.  This 
cyclic  pitch  value  is  considered  the  initial  pitch  The  thrust  moment  is 

calculated  using  the  M-flle  micalc.m.  An  initial  partial  derivative  of  cyclic  pitch  with 
respect  to  thrust  moment  is  determined  by  collectively  increasing  cyclic  pitch  a  very 
small  amount  and  determining  the  corresponding  thrust  moment.  Giving  the  new 
cyclic  pitch  and  corresponding  new  thrust  moment  the  labels  and  M(„e*, 
respectively,  the  derivative  is  determined  as; 
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dd  _  ^(nw)~^(int)  j27 

Next,  the  first  harmonic  terms  are  removed  from  the  original  thrust  moment 
Using  a  numerical  Fourier  analysis,  the  first  harmonic  terms  of  the  thrust  moment 
are: 


"is* "(new)  Sin  (111) 

where: 

n  =  number  of  azimuth  sectors 

=  total  thrust  moment  at  each  azimuth  angle 
^  =  azimuth  angle 

The  new  thrust  moment  is  then  calculated  as  follows: 

"(new)  =  "(int)  'M^c^OS  (l|r)  -W,^sin  (ijf ) 

The  change  in  thrust  moment  is  then  calculated  as; 

^"*"(new)  “"(int) 

and  the  new  cyclic  pitch  (d(ncw))  is  then  determined  by  the  relationship: 

®  (new)  *  ®  (int) 


(28) 


(29) 


(30) 


(31) 


(32) 


A  first  harmonic  function  is  then  enforced  on  theta  by  using  equation  18  where; 


(33) 


Using  the  new  cyclic  pitch  from  equation  32,  a  new  thrust  moment  is  calculated,  and 
subsequently  a  new  change  in  moment  is  determined.  The  previous  values  for  cyclic 
pitch  and  thrust  moment  become  the  initial  values,  and  the  process  is  repeated 
(skipping  the  small  collective  increase  used  for  the  initial  determination  of  the 
derivative)  until  there  is  only  a  small  first  harmonic  (2%  of  the  total  thrust  moment) 
left  in  the  thrust  moment. 

At  anytime  during  the  collective  or  cyclic  trimming  routines,  if  the  values 
of  the  respective  parameters  (difference  in  thrust  or  difference  in  thrust  moment) 
begin  to  diverge,  program  execution  will  halt.  Divergence  of  these  parameters 
indicate  that  the  blade  is  stalling  over  too  large  of  an  azimuth  region.  At  this  point, 
the  user  is  told  to  use  a  slower  airspeed  or  redesign  the  system. 

Once  the  cyclic  is  trimmed  and  the  collective  retrimmed,  trim.m,  using 
dmcalc.m,  then  calculates  the  differential  drag  values  for  each  azimuth  angle  (dD^) 
and  the  rotor  drag  moments  for  ea<h  a?imuth  angle  (DM^).  A  listing  for  dmcalc.m  is 
found  at  Appendix  H.  The  drag  moments  are  determined  by  summing  the  product  of 
the  differential  drag  values  and  radial  distance  from  the  center  of  the  hub  to  the  center 
of  each  differential  blade  element.  The  differential  drag  values  are  calculated  as 
follows; 
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dD^*  -|p.cdr(uj+y|)  [C^sin(<t>)  +CoCos  (<|>)  ] 


(35) 


where: 

p,  =  ambient  air  density 

c  =  rotor  blade  chord 

dr  =  blade  element  width 

Up  =  velocity  peipendicular  to  rotor  blade 

U,  =  velocity  parallel  to  rotor  blade 

Cl  =  coefficient  of  lift 

Co  =  coefficient  of  drag 

4>  =  inflow  angle 

The  contribution  of  profile  drag  in  the  tip  loss  region  is  accounted  for  in  the  drag  and 
drag  moment  equations  in  a  similar  manner  as  in  thrcalc.m.  The  equation  for  the 
differential  drag  in  the  tip  loss  region  is: 

|p.  c  (Up.,/  *  ( .  009  COS  )  (36) 

where: 

p,  =  ambient  air  density 

c  =  rotor  blade  chord 

R  =  rotor  blade  radius 

Reff  =  effective  rotor  blade  radius 

Up  Up  =  velocity  peipendicular  to  rotor  blade 
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Utiv  =  velocity  parallel  to  rotor  blade 
^  s  inflow  angle 

The  equation  for  the  drag  moment  in  the  tip  loss  region  is: 


tip 


R- 


(R-R.ff) 


(37) 


Once  the  differential  drag  values  are  determined,  the  H-force  (rotor  drag) 
and  the  location  of  the  resultant  thrust  vector  are  then  calculated.  The  H-force  is  a 
Hrst  harmonic  function  and  is  calculated  as  follows; 


H 


jbcos  (a„p) 


where: 


b  =  number  of  blades 

arpp  =  tip  path  plan  angle 

r^  —  start  of  aerodynamic  portion  of  blade 

R  =  radius  of  rotor  blade 

/3o  =  coning  angle 

The  cosine  and  sine  terms  of  the  rotor  H-force  are  determined  as: 


(38) 


«ic=||;dr,cos(*) 

"  jTi 


(39) 


•“  j-l 


(40) 


If  the  difference  between  the  initial  rotor  H-force  and  the  computed  value  is  more 


26 


tlum  20%  of  the  initial  one,  then  the  location  of  the  resultant  thrust  vector  is 
computed  and  the  rmor  trim  process  is  repeated.  If  the  difference  is  less  than  20%, 
then  the  location  of  the  resultant  thrust  vector  is  computed  and  the  rotor  trim  process 
is  repeated  only  if  the  new  location  of  the  resultant  thrust  vector  does  not  fall  within 
its  limits. 

The  location  of  the  resultant  thrust  vector  (r,-)  is  calculated  as  follows: 

rj.=  ^  (41) 

where: 

=  mean  of  the  sum  of  thrust  moments 
=  mean  of  the  sum  of  thrust 

If  the  difference  between  the  initial  tj  and  the  computed  one  is  more  than  1.5%  of  the 
original,  then  the  rotor  trim  process  is  repeated.  If  the  difference  is  less  than  1.5%, 
then  the  rotor  is  considered  trimmed. 

3.  Perf.m 

This  M-flle  calls  the  rotor  trim  routine,  calculates  various  rotor 
performance  parameters,  and  then  displays  and  saves  the  rotor  performance  data. 

After  calling  trim.m,  the  routine  calculates  rotor  power  required  and  the  resulting 
rotor  torque  using  the  relationships: 

n _ 

550 


^zotoz  (ri.p.  )  = 


(42) 


^rocor 


n 


(43) 


Where: 

DM^  ~  total  drag  moment  at  each  azimuth  angle 
n  ~  number  azimuth  sectors 
b  =  number  of  blades 
0  =  rotor  rotational  velocity 

The  routine  then  calculates  rotor  solidity  (u),  coefficient  of  thrust  over 
solidity  (Ct/u),  coefficient  of  torque  over  solidity  (Cq/u),  coefficient  of  H-force  over 
solidity  (CH/<r),  disk  loading  (D.L.),  figure  of  merit  (F.M.),  and  advancing  blade  tip 
Mach  number  (M^).  Solidity  is  calculated  as: 


o* 


be 

HR 


(44) 


where: 

b  =  number  of  blades 
c  =  blade  chord 
R  =  blade  radius 

Ct  was  previmisly  calculated  by  trim.m  and  Cq  and  Ch  are  calculated  as: 
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(45) 


^ZOtOT 


^disk  P«  ^tip  ^ 


c„-- 


H, 


zocoz 


^diak  P«  ^tip 


where: 

Qfotor  =  rotor  torque 
Ajii  =  area  of  rotor  disk 
p,  =  ambient  air  density 
=  blade  tip  ^)eed 
Hfotor  =  rotor  H-force 

Ct,  Cq,  and  Ch  are  then  divided  by  <r  to  get  the  CJa  relationships, 
figure  of  merit  are  determined  by: 


D.L. 


*  ZOCOZ 


F.M. 


T 

D.L. 

xocor^ 

2  Pa 

550(P,,,„) 


where: 


frotor  =  rotor  thrust 
R  =  rotor  radius 
p,  =  ambient  air  density 
Protor  =  rotor  power  (h.p.) 


(46) 


Disk  loading  and 


(47) 


(48) 


Disk  loading  and  figure  of  merit  are  set  to  zero  for  hovering  flight.  The  advancing 


tip  Mach  number  is  calculated  using  the  relationship: 


M  .  ^eipCOS(grPi>) 

49.05/?7~4r6 


(49) 


where; 


V,ip  =  rotor  blade  tip  speed 
a-iff  =  tip  path  plane  angle 
V«  =  airspeed 

T,  =  ambient  temperature  (°F) 

Once  the  rotor  performance  parameters  are  calculated,  perf.m  then  asks  the 
user  if  they  want  the  rotor  performance  parameters  displayed  on  the  screen.  If  the 
answer  is  yes,  the  data  is  presented  as  two  screens  of  input  data  and  two  of  calculated 
data  with  a  user  controlled  pause  between  screens. 

Finally,  perf.m  saves  the  rotor  performance  data  to  disk  as  two  separate 
files.  The  first  file  contains  the  single  valued  input  and  calculated  parameters  and  is 
given  the  file  name  the  user  defined  in  janrad.m  with  a  ".prf  extension.  This  file  is  a 
MS-DOS*  text  file.  The  file  is  formed  using  the  MATLAB*  DIARY  function  which 
records  a  series  of  MATLAB*  FPRINTF  commands  and  then  stores  them  to  disk  as  a 
file  named  diary.  A  call  to  the  operating  system  is  then  made  renaming  diary  to  the 
appropriate  ".prf'  file  name.  The  second  file  contains  the  vector  and  matrix  valued 
data  and  is  given  the  user  defmed  file  name  concatenated  with  an  "_p"  and  a  ".mat" 
extension.  This  file  is  a  binary  file  which  the  user  can  access  using  the  MATLAB* 
LOAD  command.  The  file  is  created  using  the  MATLAB*  SAVE  command. 
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4.  Stab.m 


The  M-fUe,  stab.m,  performs  the  helicopter  stability  and  control 
calculations,  formats  the  stability  and  control  output,  and  saves  that  output  to  disk. 
This  module  uses  the  following  M-files:  1)  cbodygrp.m\  2)  cmdbwplc.m-,  3) 
cmdbwpUi.m-,  4)  cmrgrp.m\  5)  cruise.m;  6)  ctrgrp.m\  7)  dctplots.m\  8)  dctplott.m-,  9) 
derivl.m;  10)  hmrgrp.m\  11)  hover.m;  12)  htrgrp.m\  and  13)  stabout.m.  For  details 
on  stab.m  and  the  various  routines,  See  reference  2. 
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m.  USER  INSTRUCTIONS 

A.  GENERAL 

Program  installation,  guidance  for  input  requirements,  and  rotor  performance 
output  are  discussed  in  the  subsequent  paragraphs.  These  instructions  assume  that  the 
reader  has  a  working  knowledge  of  the  MS-DOS*  operating  system  and  of  386- 
MATLAB*  version  3.5. 

B.  INSTALLATION 
To  install  JANRAD: 

1 .  run  MATLAB* 

2.  ensure  you  are  in  your  working  directory 

3.  from  the  MATLAB*  command  prompt,  type  »ldrive:  (A:  or  B:)  then 
type  "INSTALL.  The  M-file,  install.m,  will  create  a  subdirectory  named  JANRAD. 

It  will  then  copy  the  JANRAD  files  to  that  directory.  Once  the  files  are  copied,  the 
current  directory  is  changed  to  JANRAD  and  the  program  is  ready  to  run. 

C.  PROGRAM  EXECUTION 
1.  Input 

To  run  JANRAD,  at  the  MATLAB*  command  prompt,  type  »JANRAD. 

At  the  initial  screen,  enter  a  1  if  you  want  to  edit  a  previously  stored  data  file, 
otherwise  enter  a  2.  If  you  chose  to  edit  a  file,  you  will  be  prompted  to  enter  the  file 
name  without  the  extension.  If  JANRAD  cannot  find  the  file  you  entered  in  the 
current  directory,  it  will  allow  you  to  try  again.  Check  your  spelling  or  insure  you 
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are  in  the  JANRAD  directory.  Once  the  file  is  loaded,  the  edit  menu,  as  seen  in 
Figure  1 ,  will  be  displayed. 


***  EDIT  MENU  *** 

1. 

pressure  altitude 

2. 

ten^erature 

3  . 

airspeed 

4  . 

gross  weight 

5. 

number  of  blades 

6. 

blade  radius 

7  . 

blade  chord 

8  . 

hinge  offset 

9. 

blade  grip  length 

10. 

blade  twist 

11. 

blade  weight 

12 . 

#  blade  elements 

13  . 

rotational  velocity 

14  . 

#  azimuth  sectors 

15. 

lift  curve  slope 

16  . 

airfoil 

17  . 

collective  pitch 

18. 

flat  plate  area 

19  . 

vert  projected  area 

20. 

wing  area 

21. 

wing  spem 

22  . 

wing  CL 

23. 

wing  CDo 

24  . 

wing  efficiency  factor 

25. 

horizontal  tail  area 

26. 

horizontal  tail  span 

27. 

horizontal  tail  CL 

28. 

horizontal  tail  CDo 

29  . 

vertical  tail  area 

30. 

vertical  tail  span 

31. 

vertical  tail  CL 

32. 

vertical  tail  CDo 

33  . 

auxiliary  thrust 

0. 

NO  CHANGES 

Input  the  parameter  to  chamge 

Figure  1.  Edit  Menu 


Select  the  number  of  the  parameter  you  wish  to  change.  The  current  value  for  that 
parameter  will  be  displayed  followed  by  a  prompt  to  enter  the  new  value.  If  you 
choose  not  to  change  the  value,  simply  press  <  Enter  >  and  the  previous  value  will  be 
maintained.  Once  a  value  is  entered,  or  <  Enter  >  is  pressed,  the  editing  menu  will 
again  be  displayed.  Choose  the  number  of  the  next  parameter  to  change,  or  a  0  to 
exit  the  editing  menu.  If  you  chose  to  create  a  new  data  file,  prompts  for  individual 
parameters  will  be  displayed.  Enter  the  value  for  each  parameter  when  prompted. 

Enter  values  for  the  specified  parameters  regardless  of  editing  or  creating  a 
data  file  in  the  following  manner; 

1 .  PA  -  Enter  the  pressure  altitude  in  feet. 
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2.  temperature  -  Enter  the  temiTerature  in  degrees  fahrenheit. 

3.  airspeed  -  Enter  forward  airspeed  in  knots.  JANRAD  automatically  converts 
airspeed  to  feet/second. 

4.  GW  -  Enter  aircraft  gross  weight  in  pounds. 

5.  number  of  blades  -  Enter  the  number  of  blades  used  in  the  rotor  system. 

6.  radius  -  Enter  the  rotor  blade  radius  in  feet.  Measure  the  radius  from  the 
center  of  the  rotor  hub  to  tLe  tip  of  the  blade. 

7.  chord  -  Enter  the  average  rotor  blade  chord  in  feet. 

8.  hinge  offset  -  Enter  the  effective  hinge  offset  in  feet. 

9.  blade  grip  -  Enter  the  distance  from  the  center  of  the  rotor  hub  to  the 
beginning  of  the  aerodynamic  portion  of  the  rotor  blade  in  feet. 

10.  blade  twist  -  Enter  rotor  blade  twist  in  degrees.  You  can  enter  blade  twist 
as  either  a  positive  or  negative  value.  JANRAD  uses  the  absolute  value  of  blade 
twist. 

1 1 .  blade  weight  -  Enter  the  weight  of  a  single  rotor  blade  in  pounds. 

12.  rotor  blade  elements  -  Enter  the  number  of  rotor  blade  radial  elements  in 
feet.  JANRAD  adds  one  element  to  account  for  tip  loss. 

13.  rotational  velocity  -  Enter  the  rotor  rotational  velocity  in  radians/second. 

14.  azimuth  sectors  -  Enter  the  number  of  rotor  disk  azimuth  sectors. 

15.  lift  curve  slope  -  Enter  the  average  slope  of  the  linear  portion  of  the  rotor 
blade  airfoil’s  lift  curve,  where  the  lift  curve  is  Q  versus  alpha. 

16.  rotor  blade  airfoil  -  Enter  a  1  if  you  are  using  an  HH-02  airfoil.  Enter  a  2 


if  you  are  using  a  VR-12  airfoil.  The  HH-02  airfoil  data  is  very  close  to  the  NACA 
0012  airfoil  data.  Therefore,  enter  a  1  for  a  generic  airfoU. 

17.  collective  pitch  -  Enter  the  estimated  coUective  pitch  at  0.7  r/R  in  degrees. 
JANRAD  automatically  converts  the  value  to  radians. 

18.  flat  plate  area  -  Enter  the  aircraft  equivalent  flat  plate  area  in  square  feet. 

19.  vertical  projected  area  -  Enter  the  vertical  projected  area  in  square  feet. 

20.  wing  area  -  Enter  the  wing  area  in  square  feet.  Enter  a  0  if  the  aircraft  has 
no  wing. 

21.  wing  span  -  Enter  the  wing  span  in  feet.  If  you  included  the  segment  of  the 
fuselage  encompassed  by  the  wing  in  the  wing  area  value,  include  it  in  the  span  value 
also. 

22.  wing  Cl  -  Enter  the  expected  lift  coefficient  for  the  wing. 

23.  wing  Cdo  -  Enter  the  wing  profile  drag  coefficient. 

24.  wing  efficiency  factor  -  Enter  the  wing  efficiency  factor. 

25.  horizontal  tail  area  -  Enter  the  horizontal  tail  area  in  square  feet. 

26.  horizontal  tail  span  -  Enter  the  horizontal  tail  span  in  feet. 

27.  horizontal  tail  Cl  -  Enter  the  expected  lift  coefficient  for  the  horizontal  tail. 

28.  horizontal  tail  Cpo  -  Enter  the  horizontal  tail  profile  drag  coefficient. 

29.  vertical  tail  area  -  Enter  the  area  of  the  vertical  tail  in  square  feet. 

30.  vertical  tail  span  -  Enter  the  span  of  the  vertical  tail  in  feet. 

31.  vertical  tail  Cl  -  Enter  the  expected  lift  coefficient  for  the  vertical  tail. 

32.  vertical  tail  Cqo  -  Enter  the  vertical  tail  profile  drag  coefficient. 
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33.  auxiliary  thrust  -  Enter  the  expected  auxiliary  thrust  in  pounds. 

Once  the  edited  or  new  values  are  entered,  you  wUl  be  prompted  to  save 
the  data  under  a  designated  file  name.  The  file  name  must  use  a  combination  of  six 
letters,  numbers,  or  the  underscore.  If  more  than  six  characters  are  entered,  you  will 
be  prompted  to  try  again.  JANRAD  will  save  the  file  to  the  current  directory 
{JANRAD)  as  the  file  name  with  a  ".mat"  extension.  If  you  were  editing  a  data  file, 
press  <&iter>  if  you  want  to  save  the  file  with  the  original  file  name. 

The  execution  menu  is  displayed  next.  Enter  a  1  if  you  want  to  perform  a 
rotor  performance  analysis.  Enter  a  2  if  you  choose  to  perform  a  stability  and  control 
analysis.  Enter  a  3  to  execute  the  rotor  dynamics  analysis.  The  rotor  dynamics 
module  is  currently  under  development.  If  you  enter  a  3,  this  fact  will  be  displayed 
and  you  will  be  asked  to  enter  another  selection.  If  you  decide  to  edit  your  data 
further  or  you  want  to  create  a  new  file,  enter  a  4.  To  quit  JANRAD,  enter  a  5. 

From  any  of  the  selection  menus,  if  you  enter  a  character  other  than  what 
is  requested,  JANRAD  will  request  that  you  enter  an  appropriate  response.  If  you 
enter  a  letter  or  special  ASCII  character,  the  MATLAB*  prompt  "undefined  variable 
or  function"  will  be  displayed  along  with  a  repeat  of  the  initial  prompt.  If  you  enter  a 
number  other  than  requested,  the  initial  prompt  will  be  repeated. 

No  further  input  is  required  if  you  select  rotor  performance.  See 
Reference  2  for  details  on  the  input  required  if  you  select  stability  and  control. 

Details  for  the  rotor  dynamics  input  requirements  should  be  covered  in  a  future  thesis. 
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2.  Output 


If  rotor  performance  is  chosen  from  the  execution  menu,  you  wUl  be 
asked  if  you  want  the  rotor  performance  results  displayed  on  the  screen.  If  you 
respond  in  the  affirmative,  four  screens  will  be  displayed.  The  screens  will  look 
similar  to  those  shown  in  Figures  2  through  S.  The  data  in  these  figures  is  from  the 
example  helicopter  (Prouty’^  helicopter)  found  in  Reference  5,  page  669.  The  first 
two  screens  contain  the  input  data  and  the  second  two  contain  the  output  data. 


***  INPUT 

DATA  *** 

prouty 

Forward  velocity 

s 

110 

kts 

Temperature 

= 

59 

degs  F 

Pressure  altitude 

= 

0 

ft 

Gross  weight 

= 

20000 

lbs 

Number  of  blades 

s 

4 

Rotor  radius 

s 

30.00 

ft 

Blade  chord 

= 

2.00 

ft 

Blade  twist 

s 

-10.00 

degs 

Blade  lift  curve  slope 

= 

5.73 

Blade  weight 

as 

240.00 

lbs 

Rotational  velocity 

s; 

21.67 

rads/sec 

Blade  grip  length 

= 

2.25 

ft 

Hinge  offset 

s 

1.50 

ft 

press  any  key  to  continue... 

Figure  2.  First  Output  Screen  Onput  Data) 
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***  INPUT  DATA  CONTINUED  *** 

prouty 

Equivalent  flat  plate  area 

m 

30.00 

ft*2 

Vertical  projected  area 

M 

300.00 

ft*2 

Wing  area 

m 

0.00 

ft"2 

Wing  span 

9 

0.00 

ft 

Wing  CL 

m 

0.00 

Wing  CDo 

m 

0.0000 

Wing  efficiency  factor 

9 

0.00 

Horizontal  tail  area 

9 

18.00 

ft^2 

Horizontal  tail  span 

9 

9.00 

ft 

Horizontal  tail  CL 

9 

-0.30 

Horizontal  tail  CDo 

9 

0.0045 

Vertical  tail  area 

9 

33.00 

ft*2 

Vertical  tail  span 

9 

7.70 

ft 

Vertical  tail  CL 

9 

0.60 

Vertical  tail  CDo 

9 

0.0045 

Auxiliary  thrust 

- 

0  lbs 

press  any  key  to  continue . . . 

Figure  3.  Second  Output  Screen  (Input  Data  Continued) 


***  CALCULATED 

DATA 

*** 

prouty 

Fuselage  drag 

S 

1232 

lbs 

Rotor  drag 

9 

235 

lbs 

Wing  lift 

9 

0 

lbs 

Wing  drag 

= 

0 

lbs 

Horizontal  tail  lift 

= 

-222 

lbs 

Horizontal  tail  drag 

S 

9 

lbs 

Vertical  tail  side  force 

s 

813 

lbs 

Vertical  tail  drag 

= 

114 

lbs 

Tip  path  angle 

s 

4.48 

degs 

Rotor  coning  angle 

S 

5.69 

degs 

Location  of  mean  thrust  (r/R) 

s 

0.63 

collective  pitch  at  .7  r/R 

s 

6.98 

degs 

1st  lat  cyclic  term-Al  (deg) 

= 

2.07 

1st  long  cyclic  term-Bl  (deg) 

= 

4.32 

press  any  key  to  continue . . . 

*  Figure  4.  Third  Output  Screen  (Calculated  Data) 
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***  CALCULATED  DATA  CONTINUED  *** 
prouty 


solidity  «  0.085 

Disk  loading  «  0.00  lbs/ft'*2 

Figure  of  Merit  *  0.00 

CT/sigma  «  0.084 

CQ/sigma  «  0.0034 
CH/sigma  «  0.0010 
Tip  mach  of  the  adv.  blade  s  0.746 
Advance  ratio  *  0.285 

Rotor  thrust  required  (TPP)  *  20284  lbs 

Rotor  power  required  *  959  h.p. 

Rotor  torque  a  24347  ft-lbs 


Figure  5.  Fourth  Output  Screen  (Calculated  Data  Continued) 


Each  of  these  scieens  is  displayed  only  one  time.  Once  you  press  any  key  to  go  to 


the  next  screen,  you  cannot  return  to  the  previous  one.  A  copy  of  the  four  screens  is 
saved  to  the  current  {JANRAD)  directory  with  the  file  name  you  entered  at  the  "save 
instructions"  screen  and  concatenated  with  a  ".pif"  extension. 

The  vector  and  matrix  data  generated  by  the  rotor  performance  module  of 
JANRAD  is  saved  to  the  current  {JANRAD)  directory  with  the  file  name  you  entered 
at  the  "save  instructions"  screen  and  concatenated  with  a  "_p"  and  a  ".mat"  extension. 
The  following  vector  and  matrix  parameters  are  saved: 

1 .  r  -  a  vector  containing  the  radii  (in  feet)  to  the  center  of  each  blade  element. 
This  is  a  row  vector  containing  n  elements,  where  n  =  #  of  blade  elements  +  1 . 

2.  ^  =  a  vector  containing  the  azimuth  angles  (in  degrees)  around  the  rotor 
disk.  This  is  a  column  vector  containing  n  elements,  where  n  =  if  of  azimuth 


sectors. 

3.  Vj  =  a  vector  containing  the  induced  velocity  distribution  (in  ft/ sec)  along  the 
rotor  blade.  This  is  a  row  vector  containing  n  elements,  where  n  =  #  of  blade 
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elements  +  1, 


4.  0  =  a  vector  containing  cyclic  pitch  (in  degrees)  with  re^)ect  to  azimuth 
angle  at  0.7  r/R.  This  is  a  column  vector  containing  n  elements,  where  n  =  #  of 
azimuth  sectors. 

5.  01  =  a  vector  containing  blade  pitch  (in  degrees)  along  the  rotor  blade.  This 
is  a  row  vector  containing  n  elements,  where  n  =  #  of  blade  elements  +  1 . 

6.  a  =  a  matrix  containing  the  angle  of  attack  (in  degrees)  distribution.  This  is 
an  mxn  matrix,  where  m  =  of  azimuth  sectors  and  n  =  of  blade  elements  +1. 

7.  T^  =  a  vector  containing  the  total  thrust  (in  pounds)  at  each  azimuth  sector. 
This  is  a  column  vector  containing  n  elements,  where  n  =  #  of  azimuth  sectors. 

8.  =  a  vector  containing  the  total  thrust  moment  (in  ft-lbs)  at  each  azimuth 
sector.  This  is  a  column  vector  containing  n  elements,  where  n  =  #  of  azimuth 
sectors. 

9.  DM^  =  a  vector  containing  the  total  drag  moment  (in  ft-lbs)  at  each 
azimuth  sector.  This  is  a  column  vector  containing  n  elements,  where  n  =  ^  of 
azimuth  sectors. 

10.  dT  =  a  matrix  containing  the  differential  thrust  (in  pounds)  distribution. 
This  is  an  mxn  matrix,  where  m  =  #  of  azimuth  sectors  and  n  =  ^  of  blade 
elements  -f-l. 

1 1 .  dM  =  a  matrix  containing  the  differential  thrust  moment  (in  ft-lbs) 
distribution.  This  is  an  mxn  matrix,  where  m  =  of  azimuth  sectors  and  n  =  #  of 
blade  elements  -I- 1 . 
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12.  dD  =  a  matrix  containing  the  differential  drag  (in  pounds)  distribution. 
This  is  an  mxn  matrix,  where  m  =  #  of  azimuth  sectors  and  n  =  #  of  blade 
elements  +1. 

A  final  rcXor  performance  module  screen  is  displayed  instructing  you  on 
how  to  access  your  rotor  performance  data.  Use  the  MS-DOS*  TYPE  command  or  a 
DOS  editor  to  view  or  print  the  single  variable  data.  The  vector  and  matrix  data  can 
be  loaded  into  memory  using  the  MATLAB*  LOAD  command.  Some  example  plots, 
using  vector  and  matrix  data  from  the  Prouty  helicopter,  are  seen  below.  Figure  6  is 
a  plot  of  hover  induced  velocity  (ft/ sec)  versus  fraction  of  blade  radius  (r/R). 


02  CU  04  04  1 

i/R 

Figure  6.  Hover  Induced  Velocity  vs.  Blade  Radius 
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A  plot  of  total  thrust  moment  (ft-lbs)  versus  azimuth  angle  (degrees)  for  tltt  Prouty 
helic(^)ter  at  125  knots  is  shown  in  Figure  7. 


xlh* 


Figure  7.  Total  Thrust  Moment  vs.  Azimuth  Angle 

A  three  dimensional  plot  showing  the  differential  thrust  (lbs)  distribution  over  blade 
radius  versus  azimuth  angle  is  at  Figure  8.  Figure  9  shows  the  angle  of  attack 
(degrees)  distribution.  The  data  for  both  plots  is  from  the  Prouty  helicopter  at  125 
knots. 
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Figure  8.  Differential  Thrust  Distribution 


Figure  9.  Angie  of  Attack  Distribution 
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The  plcMs  in  Figures  6  through  9  were  generated  from  the  data  contained  in  a  data  file 
named  prouty _p.mat.  The  data  was  plotted  using  MATLAB* 
for  WINDOWS*,  version  4,0. 

The  preceding  paragraphs  discussed  the  details  for  the  rotor  performance 
output,  the  details  for  stability  and  control  output  are  contained  in  Reference  2.  The 
details  for  the  blade  dynamics  module  will  be  found  in  a  future  thesis. 

D.  PROGRAM  UPDATES 

JANRAD  is  written  entirely  in  the  386-MATLAB*  version  3.5  programming 
language.  This  version  of  JANRAD  is  designated  version  1 .0.  The  text  in  this 
document  corresponds  to  the  software  loaded  on  a  master  3.5”  diskette  labeled 
JANRAD  version  1.0,  September  1993.  Subsequent  changes  to  the  software  should  be 
given  a  new  numerical  designation  and  loaded  on  a  separate  diskette.  Each  new 
version  of  MATLAB*  utilizes  a  slightly  different  command  set.  As  a  result,  slight 
modifications  may  be  required  to  run  JANRAD  version  1.0  on  later  versions  of 
MATLAB*.  Additionally,  JANRAD  version  1.0  was  developed  for  use  on  computers 
using  the  DOS  type  operating  system.  JANRAD  makes  some  calls  to  the  operating 
system  using  the  MATLAB*  "!"  function.  For  this  reason  and  because  actual 
MATLAB*  commands  differ  slightly,  modifications  would  be  required  to  run 
JANRAD  version  1.0  on  an  operating  system  other  than  (i.e.,  UNIX* ).  The  ".mat" 
data  files  created  by  JANRAD  version  1 .0  can  be  transferred  to  other  operating 
systems  using  a  file  transfer  program.  Those  files  could  then  be  loaded  by  the 
particular  q)erating  system’s  version  of  MATLAB*  and  used  for  plotting. 


IV.  RECOMMENDATIONS 

Improvements  to  this  program  code  can  be  made  in  several  areas.  First,  the 
blade  dynamics  module  should  be  developed  as  a  follow-on  thesis  or  series  of  AE 
3402  (Helicopter  Dynamics)  class  projects.  As  another  thesis  project,  the  program 
code  should  be  rewritten  to  woilc  with  the  UNIX*  version  of  MATLAB*.  The  UNIX* 
version  runs  faster  and  has  enhanced  graphic  capabilities.  The  program  code  should 
also  be  updated  to  include: 

1 .  The  capability  to  offset  the  center  of  gravity  from  the  center  of  the  rotor  hub 
during  the  rotor  trimming  process. 

2.  A  method  of  incorporating  Mach  number  dependency  in  the  rotor  blade 
airfoil  data.  This  may  require  the  use  of  "look-up"  tables  versus  the  current  use  of 
parametric  equations. 

3.  A  more  accurate  method  for  determining  rotor  induced  velocity  in  forward 
flight.  This  would  be  especially  beneficial  for  airspeeds  less  than  100  knots. 

4.  Changes  to  the  code  to  account  for  the  effects  of  rotor  downwash  on  the 
wing  angle  of  attack. 

5.  Changes  to  the  code  to  incorporate  a  non-rectangular,  non-uniform,  flexible 
rotor  blade. 

These  updates  to  the  code  would  expand  the  design  students  analysis  capabilities  as 
well  as  improve  the  overall  accuracy  of  the  program.  This  would  be  especially  true 
when  utilizing  the  enhanced  capabilities  of  the  UNIX*  version  of  MATLAB*. 
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It  is  also  recommended  that  this  version  of  JANRAD  (version  1 .0)  be  kq)t  on  a 
master  disk^e  and  that  future  updates  be  ^ored  on  a  separate  disk.  This  will  ensure 
that  a  working  copy  of  the  program  is  always  available.  Thesis  level  updates  to 
JANRAD  should,  at  a  minimum,  refer  to  this  document.  Updates  that  ate  not  large 
enough  to  justify  thesis  work  should  be  documented  as  attachments  to  this  document. 
It  is  recommended  that  all  update  documentation  be  maintained  (chronologically)  in  a 
binder  along  with  this  document. 
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APPENDIX  At 


JANAAD.M 


%  Joint  Army  Navy  Rotorcraft  Analysis  emd  Design 
%  (JANRAD) 

%  Version  1 . 0 

%  September  1993 

%  MAJ  Bob  Nicholson 

%  MAJ  Walter  Wirth 


%  This  program  is  an  interactive  preliminary  design  tool 
%  developed  to  aid  the  design  student  in  determination  of 
%  initial  rotorcraft  configurations  and  in  the  calculation 
%  of  perfomumce,  stability  and  control,  emd  other  parameters. 

%  The  program  will  work  for  conventional  or  confound  rotorcraft. 

%  It  will  provide  accurate  data  for  airspeeds  less  than  10 
%  knots  emd  greater  than  or  equal  to  50  knots. 

%  VarieUSle  List  for  janrad.m,  trim.m,  thrcalc.m,  tmcalc.m,  dmcalc.m, 
%  hh02clcd.m,  vrl2clcd.m,  and  perf .m 


% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

\ 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

\ 

% 

% 


a 

Adisk 

Afh 

Afv 

airfoil 

alpha 

alphaT 

b 

B 

betao 

betat 

bhoriz 

bvert 

bwing 

cblade 

CD 

CDohori z 

CDovert 

CDowing 

CDhoriz 

CDvert 

CDwing 

CH 

CH_sig 

CL 

CLhoriz 

CLvert 

CLwing 

CQ 

CQ_sig 

CT 

CT  sig 

dD~ 

ddD 

ddDM 

ddM 

ddT 

delM 


lift  curve  slope  of  rotor  system  airfoil 
area  of  rotor  disk 

fuselage  equivalent  flat  plate  drag  area 

vertical  projected  area  (fuselage  area  under  disk) 

rotor  system  airfoil  type  {HH02/VR12) 

angle  of  attack,  rotor  blade  radial  segment 

rotor  tip  path  pl2me  2mgle 

number  of  rotor  blades 

tip  loss  parameter 

rotor  coning  angle 

geometric  angle,  rotor  blade  radial  segment 

span,  horizontal  tail 

span,  vertical  tail 

span ,  wing 

chord,  rotor  blade 

drag  coefficient,  rotor  blade  radial  segment 

profile  drag  coefficient,  horizontal  tail 

profile  drag  coefficient,  vertical  tail 

profile  drag  coefficient,  wing 

drag  coefficient,  horizontal  tail 

drag  coefficient,  vertical  tail 

drag  coefficient,  wing 

rotor  H- force  coefficient 

CH/solidity 

lift  coefficient,  rotor  blade  radial  segment 

lift  coefficient,  horizontal  tail 

lift  coefficient,  vertical  tail 

lift  coefficient,  wing 

rotor  torque  coefficient 

CQ/solidity 

rotor  thrust  coefficient 
CT/solidity 

differential  drag,  rotor  blade  radial  segment 
differential  drag,  rotor  blade  tip 
differential  drag  mcxnent,  rotor  blade  tip 
differential  thrust  moment,  rotor  blade  tip 
differential  thrust,  rotor  blade  tip 
chamge  in  total  thrust  moment 
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% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

\ 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

» 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 


Dftotal 

Dfuse 

DL 

dM 

DM^si 

dr 

Drotor 

dT 

Dhoriz 

dthetadM 

Overt 

Owing 

e 

ewing 

filename 

FM 

grip 

GW 

Hrotor 

lamdaT 

Lf total 

Lhoriz 

Overt 

Owing 

MIC 

Mis 

Machtip 

mblade 

l^si 

mu 

naz 

nbe 

omega 

PA 

phi 

phi tip 
Protor 
psi 

q 

Qrotor 

r 

R 

Rbar 

RbarT 

Reff 

rho 

rT 

solidity 

Shoriz 

Svert 

Swing 

T 

Taux 

temp 

theta 

thetalc 

thetalt 

thetao 

T^si 

twist 

op 

l^tip 


resultant  of  fuselage  drag  and  aux  thrust 
total  drag  generated  by  non- rotor  bodies 
disk  loading 

differential  thrust  moment,  rotor  blade  radial  seg 
total  blade  drag  moment  at  specific  azimuth  angle 
rotor  blade  radial  segment  width 
rotor  system  drag 

differential  thrust,  rotor  blade  radial  segment 
drag,  horizontal  tail 

change  in  cyclic  pitch  with  change  in  thrust  moment 
drag,  vertical  tail 
drag,  wing 

effective  hinge  offset 
wing  efficiency  factor 
name  of  input  file 
figure  of  merit 

length  of  inner  non -aerodynamic  portion  of  blade 
aircraft  gross  weight 
rotor  H- force 

forward  flight  induced  velocity  parameter 
total  lift  generated  by  non-rotor  bodies 
lift,  horizontal  tail 
lift,  vertical  tail 
lift,  wing 

first  harmonic  (cosine)  thrust  moment  coefficient 
first  harmonic  (sine)  thrust  moment  coefficient 
Mach  nvimber  at  rotor  blade  tip 
mass  of  rotor  blade 

total  blade  thrust  moment  at  specific  azimuth  angle 

advance  ratio 

number  of  azimuth  sectors 

number  of  blade  elements 

rotor  rotational  velocity 

pressure  altitude 

inflow  angle,  rotor  blade  radial  segment 

inflow  angle,  rotor  blade  tip 

power  required  by  rotor 

azimuth  euigle 

dynamic  pressure 

rotor  torque 

radius,  rotor  blade  radial  segment 
rotor  blade  radius 
Reff -e 
rT*Rbar 

effective  rotor  blade  radius  (tip  loss) 

ambient  air  density 

location  of  resultcuit  thrust  vector 

solidity 

area,  horizontal  tail 
area,  vertical  tail 
area,  wing 
rotor  thrust 
auxiliary  thrust 
ambient  air  temperature 
cyclic  pitch 

first  harmonic  (cosine)  of  cyclic  pitch 
first  harmonic  (sine)  of  cyclic  pitch 
collective  pitch  at  .7  r/R 

total  blade  thrust  at  specific  azimuth  angle 
geometric  rotor  blade  twist 
vertical  ccxnponent  of  velocity 
vertical  con^onent  of  velocity  at  tip 
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%  nt  horizontal  cooi}onent  of  velocity 

%  Dttip  horizontal  cooqponent  of  velocity  at  tip 

%  vi  induced  velocity 

%  Vinf  forward  airspeed 

%  vtip  tip  speed 

%  wblade  weight  of  rotor  blade 


clc 

disp('  ') 
disp('  ') 

disp('  ****************************************************************/ ) 

disp('  *  *') 

disp('  *  ***  Joint  Artny/Navy  Rotorcraft  Analysis  and  Design  ***  *') 

disp('  *  (JANRAD)  *') 

disp('  *  *') 

disp ( '  *  Version  1.0  * ' ) 

disp('  *  September  1993  '*" ) 

disp('  *  *') 

disp( '  ****************************************************************/ ) 

disp ( '  ' ) 


checklxi; 

while  checkl  >  0 

%  clearing  all  variables  in  the  MATLAB  environment 

clear 
checklsl ; 

disp ( '  ' ) 

disp('  Do  you  want  to  edit  an  existing  file  or  create  a  new  one?') 

checks 1; 

while  check  >  0 

answerxinput ( '  1.  edit  existing  file  2.  create  new  file  >>  '); 


%  ***  If  editing  em  existing  file:  get  file  name,  display  edit 

%  menu,  allow  changes  to  selected  varieJsles,  and  save  under 

%  desired  file  ncune.  Loads  to  euid  saves  frcxn  current 

%  directory  as  a  .mat  file.  *** 

if  answerxsi, 

clc 

disp  ( '  ' ) 
disp ( '  ' ) 

disp('  ***  LOAD  INSTRUCTIONS  ***  ') 

disp('  ') 

disp('  A.  Input  the  name  of  the  file  to  edit.') 

disp('  B.  The  file  was  saved  in  your  previous  session') 

disp('  with  a  ".mat"  extension.') 

disp('  C.  Do  not  include  the  extension  or  (quotations.') 

disp ( '  ' ) 

disp('  ex:  dsgnl') 

flagxO; 

while  flag  <  1 

f ilenamelxinput ( '  name  of  input  file:  ','s'); 

eval ( ['flagxexist { ' ' ' , fileneunel, ' .mat' ');']) 
if  flag  <  1, 
disp('  ') 

disp('  The  file  does  not  exist,  try  again  or  <  Ctrl-C  >') 

disp('  to  exit  program.') 

end 
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end 

aval (['load  ' .filename! ] ) 
idiile  check  >  0 


***  EDIT  MENU  ***') 

1.  pressure  altitude  2.  ten9>erature ' ) 

3.  airspeed  4.  gross  weight') 

5.  number  of  blades  6.  blade  radius') 

7.  blade  chord  6.  hinge  offset') 

9.  blade  grip  length  10.  blade  twist  ') 

11.  blade  weight  12.  #  blade  elements') 

13.  rotational  velocity  14.  #  asimuth  sectors') 

15.  lift  curve  slope  16.  airfoil') 

17.  collective  pitch  18.  flatplate  area') 

19.  vert  projected  area  20.  wing  area') 

21.  wing  span  22.  wing  CL') 

23.  wing  CDo  24.  wing  efficiency  factor') 

25.  horizontal  tail  area  26.  horizontal  tail  span') 

27.  horizontal  tail  CL  28.  horizontal  tail  CDo') 

29.  vertical  tail  area  30.  vertical  tail  span') 

31.  vertical  tail  CL  32.  vertical  tail  CDo') 

33.  auxiliary  thrust') 

0.  NO  CHANGES') 

Choi ce= input ('  Input  the  parameter  to  change:  ' ) ; 

if  isempty (choice) , 
choice*0; 
end 

if  choicessl, 

disp('  ') 

PA 

tiqpsPA; 

PAz input (' Pressure  altitude  (ft):  '); 
if  iseitpty  (PA) , 

PAstnp  ; 

end 

elseif  choicess2, 
clc 

disp('  ') 
tenp 

tmpsten^ ; 

temp^input ( ' Ten^erature  (deg  F) :  '); 
if  isenqpty  (ten^) , 
tenqpstn^ ; 
end 

elseif  choicess3, 
clc 

disp('  ') 

Vinf  »Vinf /1 . 69 

tn^sVinf  ; 

Vinf s input ('Airspeed  (knots):  ')*1.69; 
if  isenqpty  (Vinf )  , 

Vinf  stnq?*l .  69  ; 
end 

elseif  choiceEs4, 
clc 


clc 
diep('  ') 
disp(' 
disp ( '  ' ) 
disp(' 
disp(' 
disp(' 
disp(' 
disp(' 
disp  (' 
disp(' 
disp ( ' 
disp ( ' 
dispT 
disp ( ' 
disp ( ' 
disp ( ' 
disp ( ' 
disp ( ' 
disp(' 
disp(' 
disp ( '  ' ) 
disp(' 


51 


di8p('  ') 

GW 

tnp-GW; 

GWaii^ut ('Aircraft  gross  weight  (lbs):  ' ) ; 
if  iseii«>ty(GW) , 

GW-tnp; 

end 

elseif  choice««5, 
clc 

disp ( '  ' ) 
b 

tnqpsb; 

bs input ('Number  of  blades:  '}; 
if  iseo^ty  (b) , 
b>tn^; 
end 

elseif  choiceaa6, 
clc 

disp ( '  ' ) 

R 

tn^sR; 

R>input (' Blade  radius:  center  of  hub  to  blade  tip  (ft):  ') 
if  isen^ty  (R)  , 

Rstnp; 

end 

elseif  choicezs?, 
clc 

di8p('  ') 
cblade 
tmpscblade ; 

cbladesinput ('Blade  chord  (ft):  '); 
if  i sempty ( cblade ) , 
cblade«tnqp; 
end 

elseif  choices«8, 
clc 

di8p('  ') 
e 

tn^se; 

esinput ('Hinge  offset  (ft):  '); 
if  i8enpty(e), 
estmp; 
end 

elseif  choice=s9, 
clc 

di8p('  ') 
grip 

tin)*grip; 

gripsinput ('Non-aerodyn  inboard  portion  of  blade  (ft):  '); 
if  isen?)ty  (grip) , 
gripstoqp ; 
end 

if  grip  <  le-10, 
gripsle-10; 
end 

elseif  choice>clO, 
clc 

disp('  ') 

twist*' twi8t*57 . 3 

tinp*twitit; 

twistsinput ( 'Blade  twist  (deg):  '); 
if  isenpty (twist)  , 
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twist«ab8 (tmp) /57 . 3 ; 

twistaabs (twist) /57.3; 
snd 

slseif  choicsasll, 
clc 

dispC  ') 

^lade 

tnpawblade; 

wbladsainput ('Weight  of  aero  portion  of  blade  (lbs):  '); 
if  isenpty(PA), 
wblade«t8f)  ; 
end 

elseif  choice«al2, 
clc 

disp('  ') 
nbe 

tnpanbe; 

nbesinput (' Number  of  blade  elements:  '); 
if  isen9ty(nbe) , 
nbeatnqp; 
end 

elseif  choice3=13, 
clc 

disp('  ') 
omega 
tnqpsomega  ; 

oinegaainput(' Rotor  rotational  velocity  (rad/sec):  '); 
if  isen^ty  (omega)  , 
omegas tmp ; 
end 

elseif  choicessl4, 

diap('  ') 
naz 

tnqpsnaz  ; 

nazsinput ('Number  of  azimuth  sectors:  '); 
if  isenqpty  (naz)  , 
nazsto^; 
end 

elseif  choicesslS, 
clc 

disp ( '  ' ) 
a 

tn^sa; 

asinput ('Average  lift  curve  slope  (CL  vs  alpha);  '); 
if  isenqpty(a), 
astirp; 
end 

elseif  choice*«16, 
clc 

disp('  ') 
airfoil 
tnqpsairf  oil  ; 
f lagsl ; 

while  flag  >  0 

airfoilsinput ('Airfoil  1.  HH-02  2.  VR-12  >>  ' ) ; 

if  isen^ty (airfoil) , 
airfoilstn^; 
end 

if  airfoilssl, 
f lag=0 ; 
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elseif  airfoil>«2, 
flag-O; 

elae 

disp('  ') 

di8p('  ***  Enter  a  1  or  2  **•*) 
end 

end 

elseif  choice»l7, 
clc 

disp ( '  ' ) 

thetao«thetao*S7 . 3 
tnpathecao; 

thetaoainput(' Collective  pitch  at  .7  r/R  (deg):  ')/57.3; 
if  isenpty (thetao) , 
thetao«tinp/57 . 3  ; 
end 

elseif  choicesais, 
clc 

disp ( '  ' ) 

Afh 

tnpaAfh; 

Afha input ( 'Aircraft  equivalent  flatplate  area  (ft'^2)  :  '); 
if  isenpty  (Afh)  , 

Afhatnp; 

end 

elseif  choiceaald, 

elc 

disp('  ') 

Afv 

tnp*Afv; 

Afvainput ('Vertical  projected  area  (ft'‘2)  :  '); 
if  isenpty  (Afv) , 

Afvatnp; 

end 

elseif  choice««20, 
clc 

disp ( '  ' ) 

Swing 
tnpaSwing  ; 

Swings  input  (' Wing  area  (ft'*2)  :  '); 
if  isenpty (Swing)  , 

Swingstnp  ; 
end 

if  Swing  <  le-10, 

Swingsle-lO; 

end 

elseif  choicess21, 
clc 

disp('  ') 

bwing 

tnpsbwing; 

bwingsinput ('Wing  span  (ft):  '); 
if  isenpty (bwing)  , 
bwingstnp ; 
end 

if  bwing  <  le-10, 
bwingsle-10; 
end 

elseif  choicess22, 
clc 

disp ( '  ' ) 

CLwing 
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r 


tofiaCLwing; 

CLwing>ii^t ( ' Bxp«cted  CL  for  the  wing: 
if  ieenpty (CLwing) , 

CLwingstnp; 

end 

eleeif  choicessSS, 
clc 

diap('  ') 

CDowing 

tinp>CDowing; 

CDowing« input (' Wing  profile  drag  coef  (CDo) :  ' ) ; 
if  iempty  (CDowing) , 

CDowing-tng); 

end 

eleeif  choice««24, 
clc 

di^C  ') 
ewing 
tnpaewing ; 

ewing> input ('Wing  efficiency  factor  (e) :  '); 
if  isenpty (ewing) , 
ewingatop ; 
end 

if  ewing  <  le-10, 
ewingsle-lO; 
end 

eleeif  choice«B25, 
clc 

diBp('  ') 

Shoriz 

tnpsShoriz; 

Shorix=input(' Horizontal  tail  area  (ft*2) :  '); 
if  ieenqpty  (Shoriz)  , 

Shorizstn^  ; 
end 

if  Shoriz  <  le-10, 

Shorizsle-lO; 

end 

eleeif  choice>x26, 
clc 

diep('  ') 

bhoriz 

tinpsbhoriz; 

bhorizsinput ('Horizontal  tail  epan  (ft):  ' ) ; 
if  ieenpty (bhoriz) , 
bhorizstn^; 
end 

if  bhoriz  <  le*10, 
bhorizsle-lO; 
end 

eleeif  choicesx27, 
clc 

diep ( '  ' ) 

Clihoriz 
tnpxCLhoriz ; 

CLhoriz*input('E3qpected  CL  for  the  horizontal  tail:  '); 
if  ieenpty (CLhoriz) , 

CLhorizstnp; 

end 

eleeif  choicess28, 
clc 

diep ( '  ' ) 
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CDohoris 
tiqpaCDohoris  ; 

CDohoris>input ('Horizontal  tail  profile  drag  coef  (CDo) :  '); 
if  iseaqpty  (CDohoriz) , 

CDdboriz>tinp; 

end 

elseif  choice««29, 
clc 

diap  ( '  ' ) 

Svert 
tn^-Svert  ; 

Svert Bit^ut ('Vertical  tail  area  (ft'*‘2)  :  *); 
if  iBeiif>ty  (Svert) , 

SvertBtoqp; 

end 

if  Svert  <  le-10, 

SvertBle-10; 

end 

elseif  choiceBBlO, 
clc 

dispC  ') 
bvert 
tmpsbvert ; 

bvertainput ('Vertical  tail  span  (ft):  '); 
if  isenqpty  (bvert) , 
bvertatnp; 
end 

if  bvert  <  le-10, 
bvertale-lO ; 
end 

elseif  choiceaa31, 
clc 

dispC  ') 

CLvert 
tmpaCLvert ; 

CLvertainput ('Expected  CL  for  the  vertical  tail:  ' ) ; 
if  isen^ty (CLvert) , 

CLvertatop; 

end 

elseif  choiceaa32, 
clc 

dispC  ') 

CDovert 
tnpaCDovert ; 

CDovertsinput ('Vertical  tail  profile  drag  coef  (CDo):  ' ) ; 
if  isenpty (CDovert) , 

CDovertatnqp  ; 
end 

elseif  choice=s33, 
clc 

disp ( '  ' ) 

Taux 

tiii>aTaux; 

TaiiXsinput  ('Auxiliary  thrust  (lbs);  '); 
if  isenpty (Taux) , 

Tauxatnf); 

end 

elseif  choiceasO, 
clc 

dispC  ') 
dispC  ') 
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***  SAVE  INSTRDCTIONS  ***') 


diep(' 
diapi'  ') 

disp('  A.  Save  the  new  data  to  a  specified  file  name.') 

disp('  B.  Do  not  use  an  extension  or  quotations.') 

disp('  C.  Use  letter /numlDer  combinations  of  6  characters  or  less.') 

dispi'  D.  The  file  will  be  saved  with  a  ".mat”  extension.') 

disp ( '  ' ) 

disp ( '  ex :  dsgn_2 ' ) 

disp { '  ' )  ~ 

disp('  B.  If  you  made  no  changes,  press  <  Enter  >,  the  file  will') 

disp ( '  be  saved  with  the  original  name . ' ) 

f legal; 

while  flag  >  0 

f ilenameOaf ilenamel ; 

filenamelainput ('  save  file  as:  ','s'); 

if  isempty (f ilenamel) , 
filenamelaf ilenameO; 

end 

clear  filenameO 
if  length (f ilenamel)  >  €, 
disp ( '  ' ) 

disp('  use  6  characters  or  less') 

f lagai ; 

else 

flagaO; 

end 

end 

eval ( [ ' save  ' , f ilenamel] ) 
checlcaO; 
else 

disp ( '  ' ) 

disp{'  Enter  a  displayed  number  ...press  any  )cey  to  continue') 

pause 
end 
end 


%  ***  If  creating  a  new  file:  get  input  for  required  variaUsles 

%  and  save  under  desired  file  name.  Saves  to  current 

%  directory  as  a  .mat  file.  *** 

elseif  euiswer==2, 
clc 

PAainput (' Pressure  altitude  (ft):  '); 
while  isempty (PA) , 
disp('  ') 

disp r You  must  enter  a  numerical  value') 

PAainput (' Pressure  altitude  (ft):  '); 
end 

tenqpainput ('Ten^erature  (deg  F)  :  '),- 
while  isenpty  (ten^)  , 
disp('  ') 

disp ('You  must  enter  a  numerical  value') 
tempaii^ut ( 'Temperature  (deg  F) :  '); 
end 

Vinf ainput  ( 'Airspeed  Omots)  :  ' )  *1 . 69 ; 
while  isenqpty  (Vinf )  , 
disp('  ') 

disp ('You  must  enter  a  numerical  value') 

Vinf  ainput  ('Airspeed  (Jcnots)  :  ')*1.69; 
end 
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GW-ii^t ('Aircraft  groas  tiaight  (lbs): 
irtiils  issnpty  (OW) , 
diBp('  ') 

disp('You  Must  sntsr  a  nusarical  valua') 

(3W-iqput ('Aircraft  gross  wai^t  (lbs):  '); 
ond 

b>ii^t ('Huabar  of  bladss:  '); 
whils  isoaipty(b}, 
disp ( '  ' ) 

dispi'You  nust  sntsr  a  msMrical  valua') 
bailout ('Runbsr  of  bladss:  '); 
end 

Raix^ut ('Blads  radius:  cantar  of  hub  to  blade  tip  (ft):  '); 
while  isenpty(R), 
disp ( '  ' ) 

disp ('You  must  enter  a  numerical  value') 

Rainput ('Blade  radius:  center  of  hub  to  blade  tip  (ft):  '); 
end 

cbladeainput (' Blade  chord  (ft):  '); 
tdiile  ieen^ty  (cblade)  , 
disp('  'T 

disp ('You  must  enter  a  numerical  value') 
cbladeai]:^>ut ( 'Blade  chord  (ft):  '); 
end 

eainput ( 'Hinge  offset  (ft);  '); 
while  iseo9ty(e), 
disp('  ') 

disp ('You  must  enter  a  numerical  value') 
eainput ( 'Hinge  offset  (ft):  '); 
end 

gripainput ( 'Non -aerodynamic  inboard  portion  of  blade  (ft):  '); 
while  isenpty (grip) , 
disp('  ') 

disp ('You  must  enter  a  numerical  value') 

gripainput ( 'Non-aerodynamic  inboard  portion  of  blade  (ft):  ') 
end 

while  grip  <  le-10, 
gripaie-10; 
end 

twistainput ('Blade  twist  (deg):  ' ) /57 . 3; , twiBt=abs (twist) ; 
while  isempty (twist) , 
di8p('  ') 

disp ('You  must  enter  a  numerical  value') 

twistainput ( 'Blade  twist  (deg):  ' ) /57 . 3 ; , twistaabs (twist) ; 
end 

wbladeainput (' Weight  of  aero  portion  of  one  blade  (lbs):  ' ) ; 
while  isen^ty (wblade) , 
disp ( '  ' ) 

dispi'You  must  enter  a  numerical  value') 

wbladeainput ( 'Weight  of  aero  portion  of  one  blade  (lbs):  '); 
end 

nbeainput ('Number  of  blade  elements:  '); 
while  isempty (nbe) , 
disp('  ') 

disp ('You  must  enter  a  numerical  value') 
nbeainput ( ' Number  of  blade  elements :  ' ) ; 
end 

cmegaainput (' Rotor  rotational  velocity  (rad/sec) :  '); 
while  isenpty (omega) , 
disp('  ') 

disp ('You  must  enter  a  numerical  value') 

omegaa input (' Rotor  rotational  velocity  (rad/ sec) :  '); 
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•nd 

nazsii^nit (' Number  of  azimuth  sectors:  '); 
while  isenpty(naz)  , 
disp('  ') 

disp('You  must  enter  a  numerical  value') 
nazsinput ( ' Number  of  azimuth  sectors :  ' ) ; 
end 

a> input (' Lift  curve  slope  of  rotor  airfoil  (CL  vs  alpha) :  '); 
while  isenqpty  (a) , 
di6p('  '} 

dispC'You  must  enter  a  numerical  value') 

a> input (' Lift  curve  slope  of  rotor  airfoil  (CL  vs  alpha):  '); 
end 

f lag.l ; 

while  flag  >  0 

airfoil>input ('Airfoil  1.  HH-02  2.  VR-12  >>  '); 

while  isenqpty  (airfoil)  , 
disp ( '  ' ) 

dispCYou  must  enter  a  numerical  value') 
airfoilzinput ( 'Airfoil  1.  HH-02  2.  VR-12  >>  ' ) ; 

end 

if  airfoil~-l  &  airfoil~s2, 
disp('  ') 

disp('  ***  Enter  a  1  or  2  ***') 
else 
flag=0; 
end 

end 

thetaozinput(' Collective  pitch  at  .7  r/R  (deg):  ')/57.3,- 
while  isempty (thetao) , 
di8p('  ') 

disp ('You  must  enter  a  numerical  value') 
thetao=input ('Collective  pitch  at  .7  r/R  (deg);  ')/57.3; 
end 

Afhs input ( 'Aircraft  equivalent  flatplate  area  (ft''2)  :  '); 
while  isenqpty  (Afh)  , 
disp('  ') 

disp ('You  must  enter  a  numerical  value') 

Afhs  input  ( 'Aircraft  equivalent  flatplate  area  (ft'*2)  :  '); 
end 

Afv^input  ( 'Vertical  projected  area  (ft'*2)  :  ' )  ; 
while  isempty (Afv) , 
disp ( '  ' ) 

disp ('You  must  enter  a  numerical  value') 

Afv=input  ( 'Vertical  projected  area  (ft‘‘2)  :  ' )  ; 
end 

Swingainput  ( '  Wing  area,  0  if  no  wing  (ft'*2)  :  ' )  ; 
while  isenpty (Swing) , 
disp ( '  ' ) 

disp ('You  must  enter  a  numerical  value') 

Swingsinput  ( '  Wing  area,  0  if  no  wing  (ft''2)  :  ' )  ,- 
end 

if  Swing-»0, 

bwings  input  ( 'Wing  spaui  (ft):  ' )  ; 
while  isenpty (bwing) , 
disp('  ') 

disp ('You  must  enter  a  numerical  value') 
bwingsinput  ('Wing  span  (ft):  ' )  ,- 
end 

if  bwing  <  le-10, 
bwing=le-10  ,- 
end 
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CLwing«iiq)ut  ( '  Bsqpected  CL  for  the  wing :  ' )  ; 
while  i8enc>ty  (CLwing) , 
disp ( '  ' ) 

disp('You  must  enter  a  nijmerical  value') 

CLwing»input  ( '  Biqpected  CL  for  the  wing :  ' )  ; 
end 

CDowing«input ( ' Wing  profile  drag  coef  (CDo) :  '); 
while  i senqpty  ( CDowing ) , 
disp ( '  ' ) 

disp ('You  must  enter  a  numerical  value') 

CDowingsinput ( ' Wing  profile  drag  coef  (CDo):  '); 
end 

ewingsii^ut (' Wing  efficiency  factor  (e) :  '); 
while  isen^ty (ewing) , 
disp('  ') 

disp ('You  must  enter  a  numerical  value') 
ewingsinput ( ' Wing  efficiency  factor  (e) :  '); 
end 

if  ewing  <  le-10, 
ewingsle-10; 
end 

else 

Swingsle-10; 

bwing«le-10; 

CLwing«0 ; 

CDowingsO; 

ewingsle-lO; 

end 

Shoriz= input  ('Horizontal  tail  area,  0  if  none  (ft'*'2)  :  '); 
while  isempty (Shoriz) , 
disp('  ') 

disp ('You  must  enter  a  numerical  value') 

Shorizsinput  ( 'Horizontal  tail  area,  0  if  none  (ft‘‘2)  :  ' )  ; 
end 

if  Shoriz-=0, 

bhorizsinput ('Horizontal  tail  span  (ft):  '); 
while  isempty (bhoriz)  , 
disp ( '  ' ) 

disp ('You  must  enter  a  numerical  value') 
bhorizsinput ( 'Horizontal  tail  spam  (ft):  ' ) ; 
end 

if  bhoriz  <  le-10, 
bhorizsle-lO; 
end 

CLhorizsinput ('Expected  CL  for  the  horizontal  tail:  ' ) ; 
while  isen^ty (CLhoriz) , 
disp('  ') 

disp ('You  must  enter  a  numerical  value') 

CLhorizsinput ('Expected  CL  for  the  horizontal  tail:  ' ) ; 
end 

CDohorizsinput ('Horizontal  tail  profile  drag  coef  (CDo):  ' ) ; 
while  isenpty (CDohoriz) , 
disp('  ') 

disp ('You  must  enter  a  numerical  value') 

CDohorizsinput ( 'Horizontal  tail  profile  drag  coef  (CDo):  ' ) ; 
end 

else 

Shorizsle-10; 

bhorizsle-10; 

CLhorizsO ; 

CDohorizsO; 

end 
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Svert>input  ('Vertical  tail  area,  0  if  none  (ft'^2)  :  '); 
while  isespty (Svert) , 
disp ( '  ' ) 

di8p('You  must  enter  a  numerical  value') 

Svertsinput ('Vertical  tail  area,  0  if  none  (ft"‘2)  :  '); 
end 

if  Svert-»0, 

bvertsinput ( 'Vertical  tail  span  (ft):  '); 
while  isempty (bvert) , 
disp('  ') 

disp ('You  must  enter  a  numerical  value') 
bvertsinput ( 'Vertical  tail  span  (ft):  '); 
end 

if  bvert  <  le-10, 
bvertsle-lO; 
end 

CLvertsinput(' Expected  CL  for  the  vertical  tail:  ' ) ; 
while  isen^Jty  (CLvert)  , 
disp ( '  ' ) 

disp ('You  must  enter  a  numerical  value') 

CLvertsinput (' Bisected  CL  for  the  vertical  tail:  '); 
end 

CDovert=input ('Vertical  tail  profile  drag  coef  (CDo) :  ' ) ; 
while  isempty (CDovert) , 
di8p('  ') 

disp ('You  must  enter  a  numerical  value') 

CDovert = input ( 'Vertical  tail  profile  drag  coef  (CDo);  ' ) ; 
end 

else 

Svert=le-10; 

bvert=le-10; 

CLvert=0; 

CDovert=0; 

end 

Tauxsinput ( 'Auxiliary  thrust  (lbs):  ' ) ; 
while  isempty (Taux) , 
di6p('  ') 

disp ('You  must  enter  a  numerical  value') 

Tauxs input ( 'Auxiliary  thrust  (lbs):  ' ) ; 
end 


A. 

B. 

C. 

D. 


clc 
disp  (' 
disp ( ' 
disp ( ' 
disp ( ' 
disp ( ' 
disp ( ' 
disp ( ' 
disp ( ' 
disp ( '  ' ) 
disp (' 
disp('  ') 
disp ( ' 

flagsl 
while  flag  >  0 
filenamelx input ( ' 

if  isenpty (f ilenamel) 
filenamelx ' dsgn_l ' ; 
end 

if  length (filenamel) 
di8p('  ') 


) 


SAVE  INSTRUCTIONS 


) 


Save  the  data  to  a  specified  file  name.') 

Do  not  use  an  extension  or  quotations . ' ) 

Use  letter/number  ccxnbinations  of  6  characters  or  less . ' ) 
The  file  will  be  saved  with  a  ".mat"  extension.') 


ex :  dsgn_l ' ) 

If  you  do  not  enter  a  name,  the  default  is  "dsgn_l"  ') 


E. 


save  file  as; 


s' ) 


>  6, 


61 


use  6  characters  or  less') 


disp(' 
flagsl; 
else 
flag-0; 
end 

end 

eval ( ['save  ' , f ilenamel] ) 
check-0 ; 
else 


disp { '  ' ) 


disp ( '  Enter 

end 

a  1  or 

2') 

end 

clc 

dxsp ( '  ' ) 
disp  ( '  ' ) 
disp ( ' 
pause ( ! ) 

***  DATA  ENTRY  COMPLETE  ***') 

check2=l; 
while  check2  >  0 

clc 

disp('  ') 
disp ( ' 
disp ( '  ' ) 

***  EXECUTION  MENU  ***') 

disp ( ' 

1. 

Rotor  Perform£uice  Analysis') 

disp ( ' 

2  . 

Sta^Dility  and  Control  Analysis') 

disp ( ' 

!. 

Rotor  Dynamics  Analysis') 

disp(' 

4  . 

Change  data ' ) 

disp ( ' 

5. 

Quit' ) 

check! =1  ; 
while  check!  >  0 

answer-input ( ' 

Enter  a  1,  2,  !,  4,  or  5 

save  ten^  checkl 

check2 

check!  filenamel 

if  answer  ==  1, 
perf 

load  temp 
check! =0 ; 

elseif  answer  ==  2, 
stcib 

load  temp 
check! =0; 

elseif  auiswer  --  !, 

%  insert  name  of  blade  dynamics  M-file  here 
load  temp 
check!=0; 
disp ( '  ' ) 

disp('  Blade  dyneunics  module  is  currently  under  development') 

disp ( '  ' ) 

disp('  press  euiy  key  to  continue...') 

pause 

elseif  answer  —  4, 
clc 

check2-0 ; 
check! =0 ; 
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elseif  answer  » 
disp('  ') 
checklsO; 
checks «0; 
Check3s0; 
else 

disp{'  ') 
end 
end 
end 
end 
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APPENDIX  Bt  TRIM.M 


%  Rotor  trim  svibroutine 
di8p('  ') 

disp('  ***  EXECUTING  ROTOR  TRIM  ROUTINE  *** 


%  ***  calculation  of  required  parameters  *** 

rho».002377* (- .000031*PA+ (- .002*ten?j+1.118) ) ; 

%  ***  first  guess  at  rotor  profile  drag  (  H  force)  *** 

if  Vinf  <  16.9, 

DrotorsO ; 
else 

DrotorsVinf* (rho/ . 002377) ; 
end 

q=0.5*rho*Vinf*2; 

Adislc=pi*R''2  ; 

Vtip=omega*R; 

Dfu8e=q*Afh; 

CDwingaCDowing-f  (CLwing*2/  (ewing*pi*  (bwing'^2 /Swing) )  )  ; 
CDhori2=CDohorii+  (CIjhoriz'‘2/  { .8*pi*  (bhori2*2/Slioriz) ) )  ; 
CDvert=CDovert+  {CIivert''2/  ( .8*pi*  (bvert'‘2/Svert) )  )  ; 
Dwing=q*CDwing*Swing ; 

Dhorizaq*CDhoriz*Shoriz ; 

0vertaq*CDvert*Svert ; 

Df totals (Dfuse+Dwing+Dhoriz+Dvert) -Taux; 
Ijwingsq*CLwing*Swing  ; 

Lhorizsq*CLhori2*Shoriz ; 

Lvertsq*CLvert*Svert ; 

Lf totalsLwing+Lhoriz ; 

alphaTsateui2  ( (Dftotal-fDrotor)  ,  (GW-Lf total) )  ; 
musVinf*cos (alphaT) /Vtip; 


%  ***  thrust  calculation  *** 

if  Vinf  <  16.9, 

T=(l+(0.3*Afv/Adis)c)  )  *GW; 
CTsT/  (Adi8)c*rho*Vtip*2)  ; 
else 

Ts (GW-Lf total) /cos (alphaT) ; 
CTsT/  (Adis)c*rho*Vtip*2)  ; 
end 


%  ***  setup  blade  radius  elements,  azimuth  elements, 

%  induced  velocity  distributions,  amd  determination 

%  of  coning  angle  euid  tip  lose  parameter  *** 

Bsl- (eqrt(2*CT) /b) ; 

ReffsB*R; 

RbarsReff -e; 

drs (Ref f -grip) /nbe; 
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r»grip : dr : Ref f - dr ; , r*r+dr/2 ; 
rTl*0.7;,%  ***  first  guess  at  rT  *** 

RbarT=rTl*Rbar ; 
inbladeswblade/32 . 17  ; 

betaoaasin  ( (T/b*RbarT-  ( .5*  (R-e)  +e)  *vfblade)  /  ( ( .  5*  (R-e)  +e)  '*2*oinega'*2*inblad 
e) )  ; 

betatstwist* (0.7- (r/R) ) ; 

psisO :360/naz :360-360/naz; ,psi>psi' /57 .3 ; 


%  ***  induced  velocity  determination  •** 

if  Vinf  <  16.9, 

A=4*pi; 

Bvx (b/2) *omega*a*cblade; 

Tv=0; 

delT=T-Tv; 

while  aU9S(delT)  >  .01*T 

thetavstwist*  (0 . 7-  (r/R)  )-»-thetao; 

C=  ( -b/2)  *cblade*omega‘‘2*r*a.*thetav; 
vi=(-Bv+sqrt(Bv^2- (4*A*C) ) )/(2*A)  ; 

dTv= (b/2) *rho* ( (omega*r) . ^2) *a. * (thetav- (vi . / (omega^r) ) ) *cblade*dr ; 
Tvssvun  (dTv)  ; 
delT=T-Tv; 
if  delT  <  0, 

thetaosthetao- 0 . 5*thetao*abs (delT/T) ; 
else 

thetaosthetao+0 . 5*thetao*2dDS (delT/T) ; 
end 

end 

else 

lamdaTsIl  -2*inu*sin  (alphaT) 

mu'*2*  (sin(alphaT)  ‘‘2+1)  -2*mu''3*sin  (alphaT)  mu*4*Bin  (alphaT)  ‘‘2-CT''2/4]  ; 
laindaTsmax(real  (roots  (lamdaT) ) )  ; 
vi=l6undaT*Vtip-Vinf*sin(alphaT) ; 
vi=vi*ones (r) ; 
end 


%  ***  first  guess  at  theta  *** 

thetalc=0 . 035* ( (0 . 0006e- 3*Vinf “2+0 . 244e-3*Vinf ) /O . 105) ; 
thetals=-0.087* ( (0 . 0006e-3*Vinf “2+0 .244e-3*Vinf ) /0.105) ; 
theta=thetao+thetalc. *cos (psi) +thetals. *sin (psi) ; 


%  ***  rotor  trimming  routine  *** 

disp ( '  ' ) 
disp ( '  ' ) 

disp('  ***  TRIMMING  COLLECTIVE  ***') 

Jc=l; 

errorOs (T* . 02 ) +1 ; 

while  2UDs(error0)  >  T*.02 
Tpsiszeros (psi) ; 
thrcalc 

error0=T- (mean (Tpsi) *b) ; 
if  errorO  <  -T*.02, 


thetao>thetao-0.35*thetao*ab8(l.S*error0/T) * (1-mu) ; 
elseif  errorO  >  T*.02, 

thetaoBthetao-fO  .35*thetao*abs  (1 .5*errorO/T)  *  (1-mu)  ; 
end 

thetasthetao-f thetalc .  *cos  (psi )  -t-thetals .  *8in  (p8i }  ; 
if  )c  >  1, 

if  ab8(error0)  >  ab8(errorl), 
clc 

di8p('  ') 
di8p('  ') 

di8p('Thi8  configuration  will  not  trim') 
di8p('Txy  a  lower  airapeed  or  a  new  deaign' ) 
di8p('  ') 

diap (' Program  execution  terminated') 
di8p('  ') 

error {'***  END  OF  PROGRAM  ***') 
end 
end 

error 1 =error 0 ; 
lc«Jc+l; 
end 


di8p('  ') 

disp('  ***  TRIMMING  CYCLIC  ***') 

tO=clock; 

)c=l; 

errorO- ( ( (T/b) *rTl* (R-grip) ) *.04)+l; 

while  errorO  >  ( (T/b) *rTl* (R-grip) )* .04 

time=etime  (cloclc,  to)  ; 
if  time  >  15, 
disp ( '  ' ) 

disp(' still  trimming  ...') 
tO=clocJc; 
end 

I^si  {:  .k)  szeros  (psi )  ; 
tmcalc 

thetas [theta  theta ( : , k) ] ; 
t^sis  [I^si  t^si  ( : ,  k)  ]  ; 

%  ***  calculation  of  initial  dthetadM  *** 

if  k  <  2, 

theta ( : , k+l ) =theta ( : , k) +0 . 25/57 . 3 ; 
t^si ( : ,k+l) szeros (psi) ; 

ksk+1 ; 

tmcalc 

ksk-1; 

dthetadMs  (theta  ( : , k+l)  -  theta ( :  ,k) )  . /  (I^si  ( :  ,k+l)  -iSpsi  ( :  ,k) )  ; 
end 

%  ***  calculation  of  M  first  harmonic  parameters  *** 

Mlcs2*sum(l^si  ( :  ,k)  .*co8  (psi) )  /naz; 

Mlss2*sum(I^si  (:  ,k)  .*8in(psi) )  /naz; 

%  ***  removal  of  first  harmonic  terms  from  J^si  *** 
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14>si  ( :  ,k+l)  ( :  ,k)  -Mlc.^cos  (psi)  -Mis.  *Bin  (psi)  ; 

dslM-%8i  ( : ,  k-fl )  -Mpsi  { : ,  k)  ; 
errorOsmax (delM) -min(delM) ; 
if  k  >  1, 

if  errorO  >  errorl, 
clc 

disp('  ') 
disp ( '  ' ) 

dispC'This  configuration  will  not  trim') 
disp ('Try  a  lower  airspeed  or  a  new  design') 
disp('  ') 

disp ('Program  execution  terminated') 
dispi'  ') 

error ('*•*  END  OF  PROGRAM  ***') 
end 
end 

errorl serrorO ; 

%  ***  calculation  of  new  theta  *** 

delM=0.5* (1-mu) *delM; 

theta ( : , k+1) stheta ( : , k) + (dthetadM. *delM) ; 

if  errorO  <=  ( (T/b) ‘rTl* (R-grip) ) * .04 , 
theta.  *3um(theta( ; ,k) .*cos(psi))/naz; 
thetals=..: •sum (theta ( :  ,k)  .*sin(psi) )  /naz; 
else 

thetalcs2*sum (theta ( : ,k+l) .*008 (psi) ) /naz; 
thetals=2  *sum ( theta ( : , k+1 ) . *sin (psi ) ) /naz ; 
end 

theta ( : , k+1) =thetao+thetalc.*coB (psi) +thetal8. *sin (psi) ; 

%  *•*  calculation  of  new  dthetadM  *** 

thetas (theta  theta ( : , k+1) ] ; 

I^sis  [l^si  I^si  ( :  ,k+l)  ]  ; 

theta ( : , k+2) stheta( ; ,k) +0.25/57.3; 

I^si ( : ,k+2) szeros (Mpsi ( : ,k+l) ) ; 

ksk+2 ; 

tmcalc 

k=k-2; 

dthetadMs  (theta (;  ,k+2)  -theta(;  ,k) )  . / (I^si  ( : , k+2)  -^^si  ( :  ,k) ) 

ksk+1 ; 

end 


disp('  ') 

di8p('  ***  ADJUSTING  COLLECTIVE  ***') 

disp ( '  ' ) 

thetastheta ( : , k) ; 
k-1; 

errorOs (T* .01) +1; 

while  abs (errorO)  >  T*.Ol 
Tpsiszeros (psi) ; 
thrcalc 

errorOsT- (meeui(Tpsi) *b) ; 
if  errorO  <  -T+.Ol, 

thetaosthetao-0 .25*thetao*abs (1 .25*errorO/T) * (1-mu) ; 
elseif  errorO  >  T*.01, 
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thetaosthetao-t-0 .25*thetao*abs  (l  .25*errorO/T)  *  (1-mu)  ; 
end 

thetasthetao-fthetalc . *co8  (psi)  -t-thetals . *8in (psi)  ; 
if  k  >  1, 

if  ab8 ( error 0)  >  ab8(errorl). 
clc 

diepC'  ') 
diepC'  ') 

diep('Thi8  configuration  will  not  trim') 
diapC'Try  a  lower  airapeed  or  a  new  design') 
disp('  ') 

di8p(' Program  execution  terminated') 
disp ( '  ' ) 

error ('***  END  OF  PROGRAM  ***') 
end 
end 

errorl»errorO ; 

Icsk-t-l  ; 
end 


%  ***  calculating  drag  moments  *** 

Dt^siszeros  (psi)  ; 
dmcalc 


%  ***  calculating  rotor  H  force  *** 

if  Vinf  <  16.9, 

Hrotor=0 ; 
dT= [dT  ddT] ; 
dD«  [dD  ddD]  ; 

*^dT» [dT  ddT] ; 
dDs  [dD  ddD]  ; 

for  i»l; length (r)+l, 

Hlc(i)  =2*suni  {dT( : ,  i)  .  *cos  (psi) )  /naz; 

His (i) z2*sum(dD( : , i) .*sin(psi) ) /naz; 
end 

Hrotor= ( ( (b*cos (alphaT) /2) * (sum (His) -sin(betao) *Bum(Hlc) ) ) +Drotor) /2 ; 
end 


%  ***  calculating  new  rT  *** 

rT2s  ( ( (mean(t^si  ( : ,  length  (I^si  (1,  : } }  -1) )  /mean  (Tpsi) )  /R)  -t-rTl)  /2  ; 


%  ***  check  rotor  drag  and  rT,  retrim  rotor  if  required  *** 

while  abs (Drotor-Hrotor)  >  0.2*Hrotor  |  abB(rTl-rT2)  >  0.015*rTl 

if  abs (Drotor-Hrotor)  >  0.2*Hrotor, 
disp ( '  ' ) 

disp('  **♦  ADJUSTING  ROTOR  DRAG  ***') 

end 

Drotor«Hrotor ; 
if  abs(rTl-rT2)  >  0.015*rTl, 
di8p('  ') 
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di»p('  ***  ADJUSTING  MEAN  THRUST  LOCATION  ***') 

end 


disp ( '  ' ) 

dispC  ***  RBTRIMMING  ROTOR  ***') 

disp('  ') 
dT«dT( : , 1 :nbe) ; 
dD*dO { : , 1 ;nbe} ; 


%  ***  recalculating  parameters  *** 

alphaT>iatan(  (Df  total -t-Drotor)  /  (GW-L£total) )  ; 
mu*Vin£*co8 (alphaT) /Vtip; 

if  Vinf  >*  16.9, 

Ta (GW-L£total) /cos (alphaT) ; 

CT.T/  (Adi8lt*rho*Vtip*2)  ; 
lamdaTa [1  •2*mu*sin (alphaT) 

mu*2* (sin (alphaT) *2+1) -2*mu*3+sin (alphaT)  mu*4*sin (alphaT) *2-CT*2/4] ; 
lamdaTaouoc  (real  (roots  (lamdaT) ) )  ; 
vi=lamdaT*Vtip-Vinf*sin (alphaT) ; 
viavi+ones (r) ; 
end 

B=l- (8qrt(2*CT) /b) ; 

Reff=B*R; 

RbarsReff -e; 

dr»  (Reff  -grip)  /nk)e ; 

r=grip:dr:Re££-dr; ,r=r+dr/2; 

RbarT=rT2*Rbar ; 

betaoaa8in( (T/b*Rba3rT- ( .5* (R-e) +e) *wblade) / ( ( .5* (R-e) +e) *2*omega*2*mblad 
e))  ; 


%  ***  trimming  collective  *** 

tOacloc)t; 

kal; 

errorOa (T+ . 02 ) +1 ; 

while  abs(errorO)  >  T*.02 
Tpsiazeros (psi) ; 
thrcalc 

errorOaT- (mean(Tpsi) *b) ; 
if  errorO  <  -T*.02, 

thetaoathetao-0 .35*thetao*ab8 (1 .5*errorO/T) * (1-mu) ; 
elseif  errorO  >  T*.02, 

thetaoathetao+0 .35*thetao*abs (1.5*errorO/T) * (l-mu) ; 
end 

thetaathetao+thetalc . *cos (psi) +thetal8 . *8in (psi) ; 
if  Ic  >  1, 

if  abs(errorO)  >  abs(errorl), 
clc 

disp('  ') 
di8p('  ') 

disp ('This  configuration  will  not  trim') 
disp ('Try  a  lower  airspeed  or  a  new  design') 
disp('  ') 

disp ( ' Program  execution  terminated' ) 
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di8p('  ') 

error ('***  END  OF  PROGRAM  ***') 
end 
end 

errorl>errorO ; 
ksk-t-1 ; 
end 


%  ***  trimning  cyclic  *** 


ksl ; 

errorOs ( ( (T/b) *rT2* (R-grip) ) * .04) +1; 

while  errorO  >  ( (T/b) *rT2* (R-grip) )* .04 

time-etime (clock, tO) ; 
if  time  >  15. 

dispC  still  trimming  — ') 
disp ( '  ' ) 
tO-clock; 
end 

t^psi  ( : ,  k)  szeros  (psi)  ; 
tmcalc 

thetas [theta  theta ( : , k) ] ; 

^^sis  [I^si  14>8i  ( : ,  k)  ]  ; 

%  ***  calculation  of  initial  dthetadM  *** 

if  k  <  2, 

theta ( : ,k+l) stheta( : , k) +0.25/57.3; 

I^si ( : , k+1 ) szeros (psi ) ; 

ksk+1 ; 

tmcalc 

k=k-l; 

dthetadMs  (theta( ;  ,k+l)  •theta( :  ,k) )  ./  (I^si  ( :  ,k+l)  -I^si  ( ;  ,k) )  ; 
end 


%  ***  calculation  of  M  first  harmonic  parameters  *** 

MlCs2*sum(I^si  ( :  ,k)  .*cos  (psi) )  /naz; 

Ml8s2*sum(l^si  ( :  ,k)  .*sin(psi) )  /naz; 

%  ***  removal  of  first  harmonic  terms  from  I^si  *** 

14>si  ( :  ,k+l)  sl^si  (;  ,k)  -Mlc.*cos  (psi)  -Ml8.*8in(p8i)  ; 
delMsI^si  ( : ,  k+1 )  ->4>si  ( : ,  k)  ; 
errorO>max(delM)  -min(delM)  ; 
if  k  >  1, 

if  errorO  >  errorl, 
clc 

disp ( '  ' ) 
disp ( '  ' ) 

disp ('This  configuration  will  not  trim') 
disp ('Try  a  lower  airspeed  or  a  new  design') 
disp('  ') 

disp ('Program  execution  terminated') 
disp ( '  ' ) 

error ('***  END  OF  PROGRAM  ***') 
end 
end 
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error l>errorO ; 

%  ***  calculation  of  new  theta  *** 

delM>0 . 5* (1 -mu) *delM; 

theta  ( : ,  h-t-l)  «theta  ( : ,  Ic)  (dthetadM.  *delM)  ; 
if  errorO  <■  ( (T/b) *rT2* (R-grip) ) * .04, 
thetalcs2*8um (theta  ( : ,  k)  .  *cos  (psi) )  /naz 
thetals>2*sum (theta ( : , k) . *sin (psi) ) /naz ; 
else 

thetalc>2*suffl (theta ( :  .k-t-l)  .*cos(psi) )  /naz; 
thetals«2'^BUffl (theta ( :  .k-»-l)  .*8in(psi) )  /naz; 
end 

theta ( :  ,k4-l)  sthetao-fthetalc.^cos (psi)-t-thetals.*sin(p8i)  ; 

%  ***  calculation  of  new  dthetadM  *** 

thetas  [theta  theta  ( : ,  k-t-l )  ]  ; 

l^sis  [I^si  ^^8i  ( :  ,k-fl)  ]  ; 

theta ( :  ,k-t-2)  stheta( :  ,k)  -t-0 .25/57.3; 

l^Bi  ( : ,  k-t-2 )  szeros  (Mpsi  (:,  k-t-l ))  ; 

ksk-t-2  ; 

tmcalc 

kxk-2; 

dthetadMs  (theta  ( : , k-t-2)  -theta  ( :  ,k) )  .  /  (t^si  ( ; ,  k-t-2)  -t^si  ( :  ,k) )  ; 
ksk-t-1; 

end 


%  ***  retrimming  collective  *** 

thetastheta ( : , k) ; 
k=l; 

errorOs (T* .01)+! ; 

while  ab8(error0)  >  T*.01 
Ipsiszeros (psi) ; 
thrcalc 

errorOsT- (mean(Tpsi) *b) ; 
if  errorO  <  -T*.01, 

thetao«thetao-0 .25*thetao*abs (1 .25*errorO/T) * (1-mu) ; 
elseif  errorO  >  T*.01, 

thetaosthetao-t-0 .25*thetao*abs  (1.25*error0/T)  *  (1-mu)  ; 
end 

theta=thetao-t-thetalc ,  *cos  (psi )  -i-thetals .  *sin  (psi )  ; 
if  k  >  1, 

if  6d}s  (errorO)  >  abs(errorl), 
clc 

disp ( '  ' ) 
disp('  ') 

disp ('This  configuration  will  not  trim') 
disp ('Try  a  lower  airspeed  or  a  new  design') 
disp ( '  ' ) 

disp ( ' Program  execution  terminated' ) 
disp('  ') 

error ('***  END  OF  PROGRAM  ***') 
end 
end 

error IserrorO ; 
ksk-fl  ; 
end 
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%  ***  recalculating  rotor  H  force  *** 

if  Vinf  <  16.9, 

HrotorsO ; 
dT»  tdT  ddT] ; 
dD«  [dD  ddD]  ; 
else 

dT«  tdT  ddT]  ; 
dD. [dD  ddD] ; 

for  i.l:  length  (r)-t-l, 

Hlc(i)  .2*8uin(dT( : ,  i)  .*cos  (psi) )  /naz; 

Hls  (i)  .2*suin(dD( : ,  i)  .*sin{psi) )  /naz; 
end 

Hrotor.  ( ( (b*co8  (alphaT)  /2)  *  (auintHla)  -aintbetao)  *8uin(Hlc) ) )  +Drotor)  /2  ; 
end 


%  ***  recalculating  rT  *** 

rTl=rT2 ; 

rT2=  ( ( (meeuid^ai  ( : ,  length  (I^si  (1,  : ) )  -1) )  /mean  (Tpsi) )  /R)  +rTl)  /2  ; 
end 


%  ***  recalculating  drag  moments  *** 

dT.dT(: ,l;nbe) ; 
dD.dD(:,l;nbe) ; 

Dt4>8i. zeros  (psi)  ; 

dmcalc 

dT=  (dT  ddT]  ; 

dD= [dD  ddD] ; 

clc 

disp ( '  ' ) 
disp ( '  ' ) 

dispC'  ***  ROTOR  TRIMMED  ***’) 

pause ( 3 ) 
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APPEIIDIX  Cx  PBKF.M 


%  ***  Performance  output  subroutine  *** 

%  ***  Clearing  all  variables  in  the  Matlab  environment  *** 

clear 

%  ***  Loading  input  and  check  parameters  *** 

load  ten^ 

eval ( [ ' load  ' , f ilenamel] ) 
clc 

disp ( '  ' ) 


disp('  ***  ROTOR  PERFORMANCE  ROUTINE  ***') 

dispr  •******») 


disp ( '  ‘ ) 
trim 

%  ***  Calculation  of  output  parameters  *** 

Protor^mean (Dt^si ) *b*omega/550 ; 

Qrotor=mean  (Dl^si)  *b; 
soliditysb^cblade/ (pi*R) ; 

CQ«Qrotor/ (Adisk*rho*Vtip*2*R) ; 

CH»Hrotor/  (Adisk*rho*Vtip'^2)  ; 

CT_8ig=CT/8olidity; 

CQ^sig-CQ/solidity; 

CH_8igaCH/solidity; 

Machtips  (Vtip*co8  (alphaT)  +Vinf )  /  (49 .05*sqrt  (tenqp-t-460) )  ; 
if  Vinf  <  16.9, 

DL=T/ (pi*R^2) ; 

FM= (T*sqrt (DL/ (2*rho) ) ) / (550*Protor) ; 
else 
DL=0  ; 

FM=0; 

end 

clc 

disp ( '  ' ) 

disp('  ***  PERFORMANCE  CALCULATIONS  COMPLETE  ***') 

disp ( '  ' ) 

disp('  Do  you  weuit  the  performance  results  displayed  on 

screen? ' ) 
f lagsl ; 

while  flag  >  0 

answersinput ( '  1.  yes  2.  no  :>>•); 

if  answer  ss  2, 
f lagsO ; 

elseif  auiswer  i. 


%  ***  output  to  screen  *** 

clc 

disp ( '  ' ) 

disp(’  ***  INPUT  DATA  ***') 
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eval ( I ' disp ( ' ' 


'  ,f ilenamel, 


disp ( '  ' ) 
fprintf ( ' 
fprintf ( ' 
fprintf ( ' 
fprintf ( ' 
fprintf ( ' 
fprintf ( ' 
fprintf ( ' 
fprintf ( ' 
if  airfoilssi, 
disp ( ' 
else 

disp ( ' 

end 

fprintf ( ' 
fprintf ( ' 
fprintf ( ' 
fprintf ( ' 
fprintf { ' 
disp ( '  ' ) 
disp  ( '  ' ) 
disp  ( '  press  2uiy 
pause 


Forward  velocity 
Ten^erature 
Pressure  altitude 
Gross  weight 
Number  of  blades 
Rotor  radius 
Blade  chord 
Blade  twist 

Blade  airfoil 

Blade  airfoil 

Blade  lift  curve  slope 
Blade  weight 
Rotational  velocity 
Blade  grip  length 
Hinge  offset 


key  to  continue . . . ' ) 


%6.0f  kts\n' ,Vinf/l .69) 

%6.0f  degs  F\n',teii?)) 

%6.0f  ft\n',PA) 

%6.0f  lbs\n',GW) 

%6.0f  \n' .b) 

%6.2f  ft\n' ,R) 

%6.2f  ft\n' , cblade) 

%6.2f  degs\n' , -l*twist*57 .3) 

HH-02' ) 

VR-12' ) 

%6.2f  \n' .a) 

%6 . 2f  lbs\n' , wblade) 

%6.2f  rads/sec\n' , omega) 
%6.2f  ft\n',grip) 

%6.2f  ft\n' ,e) 


clc 

disp ( '  ' ) 
disp ( ' 
eval { [' disp 
disp ( '  ' ) 
fprintf ( ' 
fprintf ( ' 
fprintf ( ' 
fprintf ( ' 
fprintf ( ' 
fprintf ( ' 
fprintf ( ' 
fprintf ( ' 
fprintf ( ' 
fprintf ( ' 
fprintf ( ' 
fprintf ( ' 
fprintf ( ' 
fprintf ( ' 
fprintf ( ' 
fprintf ( ' 
disp('  ') 
disp ( '  ' ) 
disp ( 'press 
pause 


***  INPUT  DATA  CONTINUED  ***') 

( ' '  ' , f ilenamel ,''')']) 


Equivalent  flat  plate  area  =  %6.2f 
Vertical  projected  area  =  %6.2f 
Wing  area  »  %6.2f 
Wing  span  s  %6.2f 
Wing  CL  =  %6.2f 
Wing  CDo  =  %6.4f 
Wing  efficiency  factor  =  %6.2f 
Horizontal  tail  area  =  %6.2f 
Horizontal  tail  span  =  %6.2f 
Horizontal  tail  CL  =  %e.2f 
Horizontal  tail  CDo  =  %6.4f 
Vertical  tail  area  =  %6.2f 
Vertical  tail  span  =  %6.2f 
Vertical  tail  CL  =  %6.2f 
Vertical  tail  CDo  =  %6.4f 
Auxiliary  thrust  =  %6.0f 


any  key  to  continue. 


ft^2\n' ,Afh) 
ft*2\n' ,Afv) 
ft'‘2\n' , Swing) 
ft\n' ,bwing) 

\n' .CLwing) 

\n' , CDowing) 

\n' , ewing) 
ft''2\n'  , Shoriz) 
ft\n' ,bhoriz) 
\n' , CLhoriz) 

\n' , CDohoriz) 
f t*2\n' , Svert) 
ft\n' ,bvert) 

\n' .CLvert) 

\n' , CDovert) 
lbs\n' ,Taux) 


clc 

disp { '  ' ) 
disp ( ' 

eval  ( ['disp(' ' 
disp('  ') 
fprintf ( ' 
fprintf ( ' 
fprintf ( ' 
fprintf ( ' 
fprintf ( ' 


***  CALCULATED  DATA  ***') 

' ,f ilenamel, '*')']) 

Fuselage  drag  *  %6.0f  lbs\n' ,Dfuse) 
Rotor  drag  =  %6.0f  lbs\n' , Hrotor) 
Wing  lift  *  %6.0f  lbs\n' ,Lwing) 
Wing  drag  =  %6.0f  lbs\n' , Dwing) 
Horizontal  tail  lift  =  %6.0f  lbs\n' , Lhoriz) 
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fprintf ( ' 
fprintf ( ' 
fprintf ( ' 
fprintf ( ' 
fprintf ( ' 


Horizontal  tail 
Vertical  tail  side 
Vertical  tail 
Tip  path 
Rotor  coning 


fprintf (' Location  of  mean  thrust 
fprintf ('  Collective  pitch  at  . 
fprintf ('  1st  lat  cyclic  term-Al 
fprintf {'1st  long  cyclic  term-Bl 
disp('  ') 
disp('  ') 

disp( 'press  any  key  to  continue., 
pause 


drag 
force 
drag 
angle 
angle 
(r/R) 
7  r/R 
(deg) 
(deg) 


.') 


%6.0f  lbs\n' , Dhoriz) 

%6.0f  lbs\n' , Lvert) 

%6.0f  lbs\n' .Overt) 

%6.2f  degs\n' ,alphaT*57 .3) 
%6.2f  degs\n' ,betao*57 .3) 
%6.2f  \n' ,rT2) 

%6.2f  degs\n' , thetao*57 . 3) 
%6.2f  \n' , thetalc*57 .3) 
%6.2f  \n' ,thetals*57.3) 


clc 

disp ( '  ' ) 

disp('  ***  CALCULATED  DATA  CONTINUED  **•') 
eval ( ['disp ( ' '  ' . filenamel, ''')']) 

disp('  ') 
fprintf ( ' 
fprintf ( ' 
fprintf ( ' 
fprintf ( ' 
fprintf ( ' 
fprintf ( ' 
fprintf ( ' 
fprintf { ' 
fprintf ( ’ 
fprintf ( ' 
fprintf ( ' 
disp('  ') 
disp ( '  ' ) 

disp  ('press  any  Icey  to  continue...') 
pause 

clc 

f lag=0 ; 
else 

disp ( '  ' ) 

disp('  Enter  a  1  or  2') 

end 
end 


solidity  (sigma)  =  %6.3f 
Dis)c  loading  s  %6.2f 
Figure  of  Meri*  =  %6.2f 
CT/sigma  =  %6.3f 
CQ/sigma  =  %6.4f 
CH/sigma  =  %6.4f 
Tip  mach  of  the  adv.  blade  =  %6.3f 
Advamce  ratio  s  %6.3f 
Rotor  thrust  required  (TPP)  =  %6 . Of 
Rotor  power  required  =  %6.0f 
Rotor  torque  =  %6.0f 


\n' , solidity) 
lbs/ft‘2\n' ,DL) 
\n' .FM) 

\n' ,CT_8ig) 

\n' ,CQ_sig) 

\n' .CH_sig) 

\n' .Machtip) 

\n' ,mu) 

lbs\n' ,T) 

h.p . \n' , Protor) 

ft-lbs\n' .Qrotor) 


%  ***  output  to  disJt  (text  file)  *** 

disp { '  ' ) 


disp ( ' 
pause ( 2 ) 
diary  on 
diary  off 
delete  diary 
diary  on 


***  Saving  data  to  dis)c 


") 


disp ( ' 

eval  ( ('disp(' ' 
disp('  ') 
fprintf ( ' 
fprintf ( ' 
fprintf { ' 
fprintf ( ' 
fprintf { ' 
fprintf ( ' 


***  RESULTS  ***') 

' .filenamel, ''')']) 

Forward  velocity  *  %6.0f  lcts\n' , Vinf /1 . 69) 
Temperature  =  %6.0f  degs  F\n'.temp) 
Pressure  altitude  s  %6.0f  ft\n',PA) 

Gross  weight  =  i6.0f  lbs\n'.GW) 

Numlser  of  blades  »  %6.0f  \n'.b) 

Rotor  radius  =  %6.2f  ft\n'.R) 
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fprintf { '  Blade  chord 
£print£('  Blade  twist 
£print£('  Blade  li£t  curve  slope 
£print£('  Blade  weight 
£print£ r  Rotational  velocity 
£print£ ( '  Blade  grip  length 
£print£ r  Hinge  o££8et 
£print£('  Equivalent  £lat  plate  area 
£print£('  Vertical  projected  area 
£print£ ( '  Wing  area 
£print£ ( '  Wing  span 
£print£ ( '  Wing  CL 
£print£ ( '  Wing  CDo 
£print£{'  Wing  e££iciency  £actor 
£print£('  Horizontal  tail  area 
£print£('  Horizontal  tail  span 
£print£{'  Horizontal  tail  CL 
£print£ ( '  Horizontal  tail  CDo 
£print£('  Vertical  tail  area 
£print£('  Vertical  tail  spam 
fprintf ( '  Vertical  tail  CL 
fprintf ( '  Vertical  tail  CDo 
fprintf ( '  Fuselage  drag 
fprintf ( ’  Rotor  drag 
fprintf ( '  Wing  lift 
fprintf ( '  Wing  drag 
fprintf  Horizontal  tail  lift 
fprintf  V  Horizontal  tail  drag 
fprintf ('  Vertical  tail  side  force 
fprintf ('  Vertical  tail  drag 
fprintf {'  Auxiliary  thrust 
fprintf  ( '  Tip  path  auigle 
fprintf ( '  Rotor  coning  angle 
fprintf (' Location  of  mean  thrust  (r/R) 
fprintf r  Collective  pitch  at  .7  r/R 
fprintf ('  1st  lat  cyclic  term-Al  (deg) 
fprintf ('1st  long  cyclic  term-Bl  (deg) 
fprintf ( '  solidity 
fprintf ( '  Disk  loading 
fprintf ( '  Figure  of  Merit 
fprintf ('  CT/sigma 
fprintf ( '  CQ/sigma 
fprintf ( '  CH/sigma 
fprintf ( '  Tip  mach  of  the  adv.  blade 
fprintf  ( '  Adveuice  ratio 
fprintf ('  Rotor  thrust  required  (TPP) 
fprintf ( '  Rotor  power  required 
fprintf ( '  Rotor  torque 
diary  off 


%6.2f  ft\n' , cblade) 

%6.2f  degs\n' , -l*twist*57 . 3) 
%6.2f  \n' ,a) 

%6.2£  lbs\n' , wblade) 

%6.2f  rad8/sec\n' .omega) 
%£.2f  ft\n',grip) 

%6.2f  ft\n' ,e) 

%6.2f  £t‘2\n' ,Afh) 

%6.2f  ft*2\n',Afv) 

%6.2f  ft‘2\n' .Swing) 

%6.2f  ft\n'.bwing) 

%6.2f  \n'.CLwing) 

%6.4f  \n',CDowing) 

%6.2f  \n'.ewing) 

%6.2f  f t'^2\n' .  Shoriz) 

%6.2f  ft\n' . bhoriz) 

%6.2f  Xn'.CLhoriz) 

%G.4f  \n' , CDohoriz) 

%6.2f  ft"2\n' .Svert) 

%6.2f  ft\n',bvert) 

%6.2f  \n',CLvert) 

%6.4f  \n'.CDovert) 

%6.0f  lbs\n' .Dfuse) 

%6.0f  lbs\n' . Hrotor) 

%6.0f  lbs\n' .Lwing) 

%6.0f  lbs\n' . Dwing) 

%6.0f  lbs\n' . Lhoriz) 

%6.0f  lbs\n' , Dhoriz) 

%6.0f  lbs\n' . Lvert) 

%6.0f  lbs\n' .Dvert) 

%6.0f  lbs\n'.Taux) 

%6.2f  degs\n' . alphaT*57 . 3) 
%6.2f  degsNn' .betao*57 .3) 
%6.2f  \n' ,rT2) 

%6.2f  degs\n' . thetao*57 . 3) 
%6.2f  \n' ,thetalc*57.3) 

%6.2f  \n' .thetals*57.3) 

%6.3f  \n' , solidity) 

%6.2f  lbs/ft“2\n' .DL) 

%G.2f  \n' .FM) 

%G.3f  \n'.(rr_sig) 

%G.4f  \n'.CQ_sig) 

%G.4f  \n'.CH_sig) 

%6.3f  \n'.Machtip) 

VG.3f  \n'.mu) 

%G.0f  lbs\n' ,T) 

%G.Of  h.p. \n' . Protor) 

%G.0f  ft-lbs\n' .Qrotor) 


clc 

disp  ( '  ' ) 
disp  ( ' 
disp  ( '  ' ) 

di sp  ( '  A . 

eval  (  ['disp('  ' 
disp ( '  ' ) 

disp('  B. 

the' ) 
disp ( ' 

eval ( [ ' f lag=exist ( 
if  flag  <  1. 


***  PERFORMANCE  OUTPUT  DATA  INSTRUCTIONS  ***') 

Single  value  data  saved  to  a  file  named:') 

" ' .filenamel, ' .prf " '')']) 

This  is  a  text  file,  use  the  TYPE  command  to  view 

file  or  use  a  text  editor  to  view/print  the  file.') 
'  ' ' . filenamel, ' .prf' ');']) 


76 


eval ( [' ! rename  diary  ' .filenamel, ' .prf ' ' ' ] ) 
else 

eval ( [' !del  ' ,£ilenamel, ' .prf' ' ') ) 
eval  ( ['  .'rename  diary  '  ,filenamel, '  .prf' ' '  ] ) 
end 


%  ***  Output  to  disk  (.mat  file  containing  matrix  varieQiles 

%  r,  psi,  vi,  theta,  betat,  alpha,  Tpsi,  I^si,  Dl^si, 

%  dT,  dM,  dD)  *** 

%  ***  Configuring  variables  for  output  *** 

theta=theta*57 . 3 ; 

betat« (betat  twist* (0.7- (Reff+ (R-Reff) /2) /R) ] *57.3; 
alphasalpham*S7 . 3 ; , alphas [alpha  zeros (psi) ] ; 

J^sisl^si  ( :  ,  length  (Mpsi  (l,  : )  )  -  l)  ; 
dMs  [dM  ddM]  ; 

P8i=psi*57 .3; 

r=[r  (R- (R-Reff) /2)  ]  ; 

vi= [vi  0] ; 


f ilename2= [f ilenamel  '_p']; 
disp ( '  ' ) 

disp('  C.  Matrix  auid  vector  data  saved  to  a  file  neuned:') 

eval ( ['disp ( ' '  " ' , filename2, ' .mat" '')']) 

disp('  ') 

disp('  D.  This  is  a  ".mat"  binary  file,  use  the  LOAD  commamd 

to' ) 


disp('  retrieve  the  data  for  plotting.') 

eval (['save  ' , filename2, '  r  psi  vi  theta  betat  alpha  Tpsi  Mpsi  DMtosi  dT 
dM  dD'] )  -t- 


disp ( '  ' ) 
disp  ( ' 
disp ( '  ' ) 
disp('  ') 
disp ( ' 
pause 


***  END  OF  PERFORMANCE  ***') 
press  any  key  to  continue  . . . ’ ) 
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APPENDIX  Dt  THRCALC.M 


%thrcalc  calculates  the  total  thrust  along  a  blade  at 
%each  azimuth  (psi)  location 

l^>zeros  (psi*r)  ; 

UtKzeros (Up) ; 
dTs zeros (Up) ; 
ddT>zeros (psi) ; 

for  isi ; length (psi) , 

Up(i.  : ) =vi . *cos (betao) +Vin£*sin (alphaT) *cos (betao) +Vinf *cos (alphaT) *sin ( 
betao) *cos (psi (i) ) ; 

Ut(i, :) =r.*omega+Vinf*cos (alphaT) *sin (psi (i) ) ; 
phi=atan2 (Up (i, : ) ,Ut (i, : ) ) ; 
alphastheta  (i)  -fbetat-phi  ; 
if  airfoilssi, 

[CL, CD] =hh02clcd (alpha) ; 
else 

(CL,CD] =vrl2clcd(alpha) ; 
end 

dT(i,  : )  =0 .5*rho*cblade*dr*  (Up  (i,  : )  .'^2+Ut  (i,  : )  .  *2)  .*  (CL.  *cos  (phi)  -CD.  *sin 
(phi)  )  ; 

Tpsi  (i)  =sum(dT(i,  : )  )  ; 

%  ***  calculations  for  tip  loss  area  *** 


UptipsVinf*8in (alphaT)  *co8  (betao)  +Vinf  *cos  (alphaT)  *sin  (betao)  •■cos  (psi  (i) 
)  ; 

Uttip= (R- (R-Reff) /2) *omega+Vinf*cos (alphaT) *sin(psi (i) ) ; 
phitip=atam2 (Uptip,Uttip) ; 

ddT(i)  =0 . 5*rho*cblade*  (R-Reff)  *  (Uptip*2+Uttip'‘2)  *  (- .  009*sin  (phitip) )  ; 
Tpsi (i) =Tpsi (i) +ddT(i) ; 
end 
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APPENDIX  E:  TMCALC.M 


%tmcalc  calculates  the  total  thrust  moment  along  a  blade 
%at  each  azimuth  (psi)  location 

cqpszeros  (psi*r)  ; 

Utszeros (Dp) ; 
dMazeros (Dp) ; 
ddMszeros (psi) ; 

for  isl ; length (psi) , 

Up{i,  : ) avi . *cos (betao) +Vinf *sin (alphaT) *cos (betao) +Vinf *cos (alphaT) *sin ( 
betao) *cos (psi (i) ) ; 

Dt(i, :) =r. *omega+Vinf*cos (alphaT) *sin (psi (i) ) ; 
phisatan2 (Dp(i, :) ,Dt(i, ;)); 
alphaatheta(i,)c)  +betat-phi; 
if  airfoilsal, 

[CL, CD] =hh02clcd (alpha) ; 
else 

[CL, CD] =vrl2clcd(alpha) ; 
end 

dM(i,  : )  =0 .5*rho*cblade.  *r*dr .  *  (Dp  (i,  : )  .'‘2+Dt  (i,  : )  .  '‘2)  .*  (CL.*cos  (phi)  -CD. 
*Bin (phi) ) ; 

I^si  (i,lc)  =sum(dM(i,  :) )  ; 

%  ***  calculations  for  tip  loss  areas  *** 


Dptip=Vinf*sin (alphaT) *co8 (betao) +Vinf*cos (alphaT) *sin (betao) *cos (psi (i) 

)  ; 

Dtcipa (R- (R-Reff ) /2) *omega+Vinf *cos (alphaT) *sin (psi (i) ) ; 
phitip=atem2 (Dptip,Dttip) ; 

ddM(i) =0 .5*rho*cblade* (R- (R-Reff) /2) * (R-Reff) * (Dptip“2+Dttip*2) * ( - . 009*s 
in (phi tip) ) ; 

I^si  =^^si  (i,]c)  -i-ddlKi)  ; 
end 


APPENDIX  Ft  DMCALC.M 


%dmcalc  calculates  the  total  drag  along  a  blade  at 
%each  azimuth  (psi)  location 

Up>zeros (psi*r) ; 

Utszeros (Up) ; 
alphams zeros (Up) ; 
dD=zeros (Up) ; 
ddDszeros (psi) ; 
ddDMszeros (psi) ; 

for  isl '.length (psi )  , 

Up (i, : ) =vi . *cos (betao) +Vinf *sin (alphaT) *cos (betao) +Vinf *cos (alphaT) *sin ( 
betao) *cos (psi (i) ) ; 

Ut (i, : ) »r.*omega+Vinf*cos (alphaT) *sin(psi (i) ) ; 
phisateui2  (Up  ( i ,  : )  ,  Ut  ( i ,  : )  )  ; 
alphastheta  ( i )  -t-betat  -phi  ; 
alpham ( i , : ) ^alpha ; 
if  airfoil==l, 

[CL.OD]  =hh02clcd (alpha)  ; 
else 

[CL, CD] =vrl2clcd (alpha) ; 
end 

dD(i, :)=0.5*rho*cblade*dr* (Up(i, :) .^2+Ut(i, :) .^2) . * (CL. *sin (phi) +CD. *cos 
(phi) ) ; 

dDM=dD(i, :) .*r; 

Dt^si  (i)  sSum(dDM)  ; 

%  ***  calculations  for  tip  loss  area  *** 


Uptip=Vinf *sin (alphaT) *cos (betao) +Vinf *cos (alphaT) *sin (betao) *cos (psi (i) 
)  ; 

Uttip= (R- (R-Reff) /2) *omega+Vinf*cos (alphaT) *sin(psi (i) ) ; 
phitip=at2ui2 (Uptip, Uttip) ; 

ddD  (i)  =0 . 5*rho*cblade*  (R-Reff)  *  (Uptip*2+Uttip'‘2)  *  ( .  009*cos  (phitip) )  ; 
ddDM(i)=ddD(i)*(R- (R-Reff) /2) ; 

Dl^si  (i)  =DI^si  (i)  ■•■ddDM(i)  ; 
end 


80 


APPENDIX  6t  HH02CLCD.M 


function  [CL, CD] Khh02clcd (alpha. Mach) 

%hh02clcd  calculates  CL  and  CD  for  an  HH-02  airfoil 

%given  angle  of  attack  (alpha)  in  radians  and  Mach  number  (Mach) 

% [CL, CD] «hh02clcd (alpha, Mach) 

CLszeros (alpha) ; 

CD-zeros (alpha) ; 
a«alpha* 1 8  0 /pi ; 

chklz(a>s20  &  a<Bl80); 

CL»CL+chkl.* (0.42541+0. 026863*a+5.5988e-4*a. ‘2-2 .1493e-5*a. *3+1. 5932e- 7* 
a.‘4-3.4659e-10*a.‘5) ; 

CD=CD+chkl.* (-0.7179+0.061213*a-5.9861e-4*a.‘2+7.3708e-6*a.‘3-6.6605e-8* 
a.‘4+1.913e-10*a.‘5) ; 

chkls (a>s-180  &  a<s-50); 

CL.CL+chkl.* (-4 .6183-0. 1923*a-3. 5SS4e-3*a. ‘2-3 .3273e-5*a. ‘3-1 .4S28e-7*a. 
‘4-2.3003e-10*a.‘5) ; 

CD=CD+chkl.* (2.7093e-2-2.1309e-2*a+2.0335e-4*a.‘2+3 .47e-7*a.‘3-3 .0586e-8 
*a.‘4-1.2584e-10*a.‘5) ; 

chkl=(a>-50  &  a<-20); 

CL*CL+chkl.* (-2. 5519-0. 22847*a-9.5667e-3*a. ‘2-1. 7051e-4*a. ‘3-1. 0909e-6*a 
.*4)  ; 

CD=CD+chkl.* (2 .7093e-2-2 .1309e-2*a+2.0335e-4*a.‘2+3.47e-7*a.‘3-3 .0586e-8 
*a.‘4-1.2584e-10*a.‘5) ; 

chkl=(a>=-20  &  a<=-10) ; 

CL=CL+chkl . * ( -  0 . 2+0 . 089*a+0 . 0034*a. ‘2) ; 

CD=CD+chkl.* (2.7093e-2-2.1309e-2*a+2.0335e-4*a.‘2+3.47e-7*a.‘3-3.0586e-8 
*a.‘4-1.2584e-10*a.‘5) ; 

chkls(a<20  &  a>-10) ; 

CL=CL+chkl  .  *  (5  .  8766e-2+l  .  3131e-l*a+2  .4742e-3*a.  ‘2-5  .  303e-4*a.  ‘3-1  .  5818e- 
5*a.‘4+1.28e-6*a.‘5) ; 
chk2sa<-4; 
chk2schk2 . *chkl ; 

CT)»CD+chk2.*(l. 3786+0. 916*a+0.21396*a. ‘2+2. 0371e-2*a. ‘3+7. 0076e-4*a. ‘4) ; 
chk2s(a>s-4  &  a<=7) ; 
chk2schk2 . *chkl ; 

CD-CD+chk2.*(9.732e-3+3.2326e-4*a+1.4392e-4*a.‘2-8.5073e-5*a.‘3+1.1826e- 
6*a.‘4+1.5271e-6*a.‘5)  ,- 
chk2>a>7; 
chk2=chk2 . *chkl ; 

CD.CD+chk2.* (1.842e-l-5.7532e-2*a+5.8043e-3*a.‘2-1.2803e-4*a.‘3) ; 
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APPENDIX  B*  VE12CLCD.M 


function  [CL.CD] Bvrl2clcd(alpha) 

%vrl2clcd  calculates  CL  and  CD  for  the  VR-12  airfoil 
%given  angle  of  attack  (alpha)  in  radians 
% [CL, CD] >vrl2clcd (alpha) 

CLszeros (alpha) ; 

CD«zeros (alpha) ; 
a>alpha*180/pi ; 

chk>(a>«20  &  a<sl80); 

CL«CL+chk.* (1.1733 -0.018879*a+l. 5762e-3*a. *2-3 .1925e-5*a. "3+2. 0949e-7*a. 
*4-4.3807e-10*a.*5) ; 

chkx(a>=-180  &  a<x-50); 

CL»CL+chk.*(-4.6183-0.1923*a-3.5554e-3*a.*2-3.3273e-5*a.*3-1.4528e-7*a.* 
4-2 .3003e-10*a.*5) ; 

chks(a>-50  &  a<-30); 

CLaCL+chk.* (-0.22114+0. 020857*a+2.8571e-4*a.  *2)  ; 
chk3(a>a-30  &  a<=-10); 

CL«CL+chk.* (-l.ll-0.12383*a-0.01515*a.*2-6.8667e-4*a.*3-le-5*a.*4) ; 
chka(a<20  &  a>-10); 

CL=CL+chk.*(0.11976+0.12341*a+5.5841e-4*a.*2-2.0652e-4*a.*3) ; 
chks(a>al7  &  a<sl80); 

CD=CD+chk.* (-0.26376+0. 017917*a+6.9927e-4*a. *2-9. 1137e-6*a. *3+2 .6277e-8* 
a. *4)  ; 

chka(a>a-180  &  a<=-10); 

CD=CD+chk.* (-0. 17486-0. 034463*a-1.0233e-4*a. *2-2 .8958e-6*a. *3-4 .6577e-8* 
a.*4-1.5557e-10*a.*5) ; 

chkx(a>-10  &  a<s0) ; 

<3)=CD+chk.*  (9  .8678e-3+3.4934e-3*a+1.4844e-3*a.*2-1.3564e-4*a.*3-1.0936e- 

5*a.*4)  ; 

chks(a>0  &  a<xl5) ; 

CHJxCD+chk.* (9.8e-3+7.0457e-4*a+5.6104e-5*a.*2-4.1151e-5*a.*3+3.8695e-6*a 

."4)  ; 

chks(a>15  &  a<17); 

CDaCD+chk.* (-1.33+1.325e-l*a-2.5e-3*a.*2) ; 
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